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ELASTO-DYNAMIC PROBLEM CONCERNING THE 
SPHERICAL CAVITY 


By A. CEMAL ERINGEN+ 
(Purdue University, Lafayette, Indiana) 
Received 15 June 1956] 
SUMMARY 
By using the Fourier transform technique the solution is given for the elasto- 
dynamic problem resulting from the application of arbitrary dynamical tractions 
over the surface of a spherical cavity in an infinite, isotropic, elastic solid. Various 
special cases are studied. Computations are carried out and plotted in the case of 


blast loading. 


1, Introduction 

ELASTIC wave propagations originating from the application of forces on 
the surface of a spherical cavity have been studied by K. Sezawa (1), K. 
Sezawa and K. Kanai (2), W. Inouye (3), H. Kawasumi and R. Yosiyama 
(4), G. Nishimura (5), and Vanék (6). The basic aim of these studies was to 
throw light on the origin of earthquake phenomena. Consequently, the 
surface traction distribution on the surface of the cavity was limited to some 
particular types. Moreover, attention was directed to the behaviour of the 
solution at great distances from the cavity. Basically the cases studied 
include the azimuthal distribution of the type P,(cos@), P,(cos@) with 
harmonic or exponential type of loading in time.{ 

A few recent papers (7—10)§ have been devoted to the study of the one- 
dimensional problem, namely, the propagation of elastic waves originating 
from the application of uniform time-dependent pressure on the surface of 
the spherical cavity. Of these, the papers by Allen and Goldsmith (7, 8) 
represent a valuable contribution to this problem on account of the exten- 
sive set of curves for the displacement, velocity, and stresses. These papers 
assume an exponentially decaying pressure in time. 

While the above studies have important application to the problem of 
waves generated by earthquakes and explosion, they exclude some impor- 
tant special cases of loading, such as blast, concentrated impact, ring 
loading, and the general case of arbitrary loading in space and time vari- 


ables. Also the behaviour of the solution in the neighbourhood of the cavity 
is intricate and requires further attention. 
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The present paper deals with the last-mentioned class of problems which 
appear to be basic theoretically and are important from the point of view of 
applications in engineering. By using the Fourier transform we give the 
general solution of the elasto-dynamic problem for a spherical cavity under 
the action of arbitrary dynamical surface tractions. The solution given in 
all of the foregoing papers may, of course, be obtained as special cases of 
the present solution. 

Various particular cases of practical importance are studied. Among 
these, the uniformly applied normal dynamic load leads to a closed form 
solution. Curiously enough, this is not the case for the analogous problem 
in two dimensions, namely, the circular hole in an infinite plate.t Computa- 
tions are carried out for the case of blast and discussed in the body of the 


paper. 


2. Formulation of the problem 





ELAST 


It is 
going | 


The equations of small motion of a homogeneous, isotropic, elastic body ' 


under no body forces are 
j 9 
(A+p)U(V.u)+pV2u = puy. (1) 
Here u is the displacement vector; A, » are the Lamé constants, p is the 
mass density, Y is the gradient operator, VY. and V? represent the divergence 


and Laplacian; ¢ is the time, and the comma a differentiation, ' 


e. Uy, = 6?u/et?. In spherical coordinates (7, 0,4) (Fig. 1), u has the com- 
ponents wu, v, w. Components of stress a. O,5+.+, 794 are related to dis- 
placement components wu, v, w by 

- 9 
» = AA+ 2p ,, 
og = AA+2ur-l(u+v,9), 


Co 


og = AA+2yr-!(u-+-v cot 6+ w 4 cosec 4), 


9 
79 = p(v,—r-lv+r- ug), ‘ 
Tr = plw,—r-lw-+ (rsin 6) u 4], 
Tog = pr-(wg—w cot 6+ v4 cosec 8), 

where A=V.u = u,4+r[2u+v9+v cot 6+ (sin @)-1w 4] (3) 


is the dilatation. 
The problem is now reduced to finding the solution of (1) which vanishes 
* = 00 and at the surface of the spherical cavity, r = a, satisfies the 





boundary conditions c,(a, 0,4, t) = R(6, 4, t), 
7,6(a, 9, , t)= st O( 6, p, t) t), (4) 
7,4(4, 9,6, t) = (8, ¢, t). 


+t These results are now being submitted for publication. 
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rhich It is further required that the waves generated at r = a will be of out- 
ewof going type. 
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is the 3, The solution 
gene We use the Fourier transform technique to solve the problem stated. 
ation, " The Fourier transform f(p) of f(t), and f(t), are related to each other by 
—_ the equations 
0 dis- x 
- . , ‘ l r =x . = 
f( p) | f(tje” dt, fO=s S(pye-" dp. (5) 
. a@7T 
Multiplying both sides of (1) by e’”‘ and integrating the result from —0oo to 
x 7} |} re i > o\= 9 9 - 
A iia (V2+-A2)a = (1—h3/h2)Y(V. a), (6) 
where h? = pp?/(A+ 2p) x*p?, h? = pp*/p = B*p'. (7) 
Taking the divergence of both sides of (6) we obtain 
(V?+h2)7.0 = 0. (8) 
Assuming that Y. i is determined from (8) a particular solution 0, of (6) is 
nishes _ s _ 
nn the u, = A; *V(V.). (9) 
A more complete solution of (6) is obtained by adding to i, the comple- 
4) Mentary solution 0, of equations 


(V2+h2)0,=0, V.a,=0. (10) 


The desired solution must satisfy the boundary conditions at r = a and 
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r = 00 and must have the character of diverging waves. The latter con. the cot 
dition is guaranteed by the Sommerfeld radiation condition, i.e. 
lim r(0,,,—th, 0,) = 0 (k = 1,2). (1 


ro 
) 


A solution of (8) may be obtained in terms of the spherical surface har. 


monies S,, of positive integral degree n in the form } Ase 
V.a, = 7"R,(r)S,, (12 
where #,,(r) is a function of r only which must satisfy the equation ere 
surtac 
Rp t+2(n+1)rk, , +h? R, = 0. (13)) ion 
The solution of this equation is given by 
R,, = r-)[AH™ (hy r)+ BH) (hy r)], (14 


where A and Bare arbitrary constants and H\) , and H®), arethe halforder The di 
Hankel functions. Of these functions only H'‘?) , satisfies the Sommerfeld: given 
radiation condition, i.e. the asymptotic form of H“., corresponds to out: spheri 
going waves since P 
H')) (hy, pe?! ~ (2/zh;, p)' exp| —i(n+ 1)7/2] exp] i(h;, p—pt)]. z- 


Hence we set B = 0. 


Also, absorbing A into S,, we write 


v.08, = ra” (hy r)S,,. (15 


The expression of 0, now follows from (9): 


Eq 
il —2 —} (1) “\ O 4 . 
B, = —AyP?V (rt AY) (4, 7)S8,,). (16) willa 
The components @,, 6,, @, of i, are therefore given by: equat 
tracti 
3 ~2f =} FF (1) Ys) 
Uy hts -h; [2 H}, 1 y(hy U )S, wr? stress 
1 aT (1) \Q " 
U7 hy *[r tA (hy 7)S,, | 9, (1%) trans: 
=, 27 .—3 FF) \ ‘ 
u 1 ‘ m hj | Eo H\, (hy U )S, | 4 cosec 0. ¢, 
We now proceed to determine a solution @, of (10). Let Pm 
a, = £,7Z., (18 
where Z,, is a vector point function having solid spherical harmonics a ag 
rectangular components, and F,, has the same form as (14) with B=! 
and h, replacing h,. The function 
a 79 
Z, = (VW,)xr a ae 
satisfies both equations (10). Here W, is a spherical surface harmonic 0, ré 
degree n. Hence a solution of (10) is of the form 706 


a, = rH) (hor)(VW,) Xr, (20 
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the components @,, 0, @, of which are: 
ii 0, 
By rH .(har)W,, | 4 cosec 8, 
iD, [rH s(hor)W, | 6. 
A second solution of (10) may be obtained by observing the fact that 
Vx, = Vx(R#, VW, xr (22) 


Hence in this equation if we replace W, by another spherical 


(21) 


u, 
satisfies (10). 


face harmonic Y, of integral degree n we obtain anathet solution of (10): 


iy = rH) (hgr)n(n+DY,, 
=r sient (hor )¥, | 0» (23) 
i = r [rH (her)¥, |g cosec 6. 
The diverging wave type of solution of (6) which vanishes at r = 00 is then 
given by the sum 0 = 0,+i0,+ 04, whose components in terms of arbitrary 
spherical surface harmonics S,,, W,, Y,, of integral degree n are 
i hp *[rtH, (hyr)S,),+7 TH y(hor)n(n+ DY, 
b= {—hy*r tH (hy r)S rt) (hor)Yn | te+ 
+[rH\) ,(hor)W,|4cosec@, (24) 
w : h; “) Ha (hy r)s rri H' 1) (hor)Y, | +4 cosec 6— 


—[r-H® ,(hr 


Hence any linear combination 


)W,, |. 


Equation (24) satisfies (6) for every n. 


will also be solution of (6). By summing over n from 0 to oo in each 
equation of (24) we can express the motion resulting from any dynamic 
traction at the surface of the cavity. This requires the computation of 


stress compe 


ments which are obtained by substituting (24) into the Fourier 


ransforms of (2). Hence 

ul A,,( pr)S, +B,,( pr)n(n+ DY, ), 

u{B,,( pr)S, +F,( pr)n(n+1)¥, +[4,pr)S, +4, (pr), | 04 

+-[rF,( pr)W,, cosec 6] 9g}, (25) 
G4 uf FE ( pr)S F ( pr)n(n+1)Y,, + (cosec?@ 6?/e¢?+- 
cot 6 @/26)| G,,( pr)S, +H,( pr)Y, |—[rF,(pr)W, cosec 6] 9}, 

7 u{[C,,(pr)S, +D,,( pry, |e+[4rB,,( pr)W, ].4 cosec 9}, 
Tr i{[C,,( pr)S, +D,,( pr)Y,, | 4 cosec 6@—[4rB,,( pr)W,, |.6}, 
yi G,,( pr)S, +H, ( pr)Y, |eosec 9} 9g 


ple? 06? +-n(n+1) 2\[rF,( pr)W, ], 
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where 
A,( pr) = r-(apr)-*{[ (Bpr)®—2(n—1)(m-+-2) JH) (apr) + 
aallles lldeaek 
B,( pr) = Se esaiaiss, (Bpr)—3H\}) ,(Bpr)], 
C,,( pr) = —r-*(apr) haa (xpr)—3H 7" (pr), 
Dal pr) = Mi Bpr)?—2n—1)(n+ 2) (Bpr + 
|+- 2BprH\"),(Bpr)—3H ,(Bpr)} 
E,( pr) = r-*(apr)-*{[ (Bpr)?—2(apr)?—2]H ,(apr)— (26 


[2(apr) Hi? (apr)—3H{) \(apr)]}, 
F,( pr) = 2r-*H_ ,(Bpr), 
G,( pr) = —2r-*(apr)?H (apr), 
H,,( pr) A (Bpr)+ 2BprH) (Bpr)|. 





ELAST' 


express 
Legenc 


| In ord 


satisfie 


This 
from ( 


. . ' 
Here primes denote differentiation with respect to the argument of the 


functions. 
The Fourier transforms of the boundary conditions (4) give: 


R= p> [A,(pa)S,+B,(pa)n(n+ Dy, ], 


© = p> {[C,(pa)S, +D,,( pa)Y,, |o+[4aB,,( pa)W, | 4 cosec 6}, 


n 


® = p> {[C,(pa)S,+D,,( pa), ] 4 cosecO—[4aB,(pa)W,] 9}. (27 


¥ 


~ 


From these series we can determine the surface harmonics S,,, Y,,, and V,. 
From (27) it is clear that a series representation of the following form is a 
convenient one: 

D D Vi | = > Ae 

RK ym O = > (Ti, pt7n,g cosec 9), 

n n 
® = > (T,, 4 cosec 0—7, 9). (28 
n 

This particular form is also seen to have the physical significance that the 
tangential tractions consist of components ¥ 7), in 6-direction and 
> 7,4 cosec @ in d-direction, which are derivable from a potential, and of 


components > 7,4 cosec 6 in 6-direction and ¥ 7, g in ¢-direction, which are | 


derivable from a stream function. Here R,,, T. and 7,, are surface spherical 
harmonics of integral degree n, of which R,, is given by: 


=" m. 
R,, sa P,(cos@)+ ¥ (a™ cosmd-+b” sin md) P?”"(cos 8), (29 


m=1 
where a” and b” are functions of p. Surface harmonics 7’, and 7,, hav! 
similar forms, with R,,, a”, and b” replaced by 7, n(n-+-1), c”", and d™ in the 
expression for 7',, and they are replaced by n(n-+- 1), « 


m and f% in the 
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expression for 7,. Given #, we can determine a7”, b” from the Fourier- 


n?> 


Legendre expansion theorem: 


’ 7 
2n+1(n—m)! f¢ 
a ; 1 | R P”"(cos 6)cos md sin 6 dédd, 
a7 8 6((” 
7 0 (30) 
2n l (n m)! a 4 5D ° . 
b . | R P*(cos #)sin m¢ sin 6 ddd. 
27 (7 m ? 
v7 U 
(26); In order to obtain 7, and 7, we use the following differential equation 
satisfied by the spherical surface harmonics 7", 
(sin @)-"(sin #7, 9) g+(sin @)-*T), 44 —n(n+-1)T,. (31) 


This of course is also valid for 7,. If we now form the left side of (31) 


m (28) we obtain 


of the S n(n 1)7' cosec 4| (© sin @) » + b4i = T (32) 
>) n(n 1)+ cosee AO 4—( ® sin 6 ) a] = +, 7 


These equations are of the same form as the first equation of (28). Hence 

using T'/n(n+-1) and 7/n(n+-1) with 7 and 7 defined in terms of © and ® 

32) in place of R in (30), we obtain the coefficients c’", d”", and e”, f™ 
respectively 

val We now substitute (27) into the first equation of (28) and into (32), the 

latter being equivalent to the remaining two equations of (28). Using 


id WW. sa d . 
31) for all surface harmonics we obtain 
m Is 
A,( pa)S B,,( payn(n+1)Y,, por, 
C,( pa)S,, +D,,( payY,, pT, (33) 
98 aB,(pa)W, Zu". 
The solution ot (33) is 
att ; 7 
' Ss (uA dD pa)R, B,,( pa)n(n- 1)T7,, |, 
l 
nd x, (uA | A pa)T, C, ( pa)k, |, (34 
el Wi %aunkb a) a? 
el where 
\ \,,( pa) 1. ( pa)D,,( pa)—n(n+-1)B,,( pa)C,,( pa). (35) 
oq. Lhis of course completes the problem, for (34) used in (24) and (25) gives 
the Fourier transforms of displacement and stress components; on multi- 
not? : , i , : 
have» Plication by (27)-1e-'’” and integrating with respect to p from —oo to « 
in tl we obtain the actual displacement and stress components. 


in th We now give some special cases which are of practical interest. 
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4. Special cases 
(i) Zero surface shear. © = ® = 0. Hence T, = 7, = 0. 
(ii) Normal traction with central symmetry and (i). In this cage 
a” — 6b” — 0 when m is odd. 
(iii) Normal traction with constant amplitude over a cap 0<6 <6 
—nr <d¢< rand (i), ie. 
5 R(p) when 0<0<0,,—ar<db<77 
R(0, 4, p) = |) —T . 
10 otherwise. 
Substituting (36) into (30) gives: 


bm = 0, a® = (2n+1)R(p)P7(cos8,)sind,, a” =0 (m0). (37 


nt 7 di 


(iv) Normal traction with constant amplitude over a zone 0, < 0 <6, 
a <¢ <7, and (i). This follows from (37) by substituting 6, for 0, in 
(37) and subtracting the result from (37). Hence 


«x @, a® (2n L1) Ri P)| P~}(cos 6, )sin 6,— P,,*(cos @,)sin 4], 


an 0 (m+~0). (38) 


n 


(v) Concentrated load at one pole and (i). In this case we have (37) with 
lim R( p)2z7a*(1—cos6,) = R,( P)s (39) 
6,—0 
k->« 


where #,( p) is the Fourier transform of the concentrated time dependent 








— 


load R,(t) acting at the pole r = a, 6 = 0. Using (37) in (39) and carrying ; 


out the indicated limit we obtain: 


2n+-1 
by = 0, = Rp), a™"—0 (m0). (40 

47a* 
(vi) Impulsive concentrated load and (i). In this case we write 
R,(t) R, d(t) where R, is constant and 8(t) is the Dirac delta function 
defined by 


s(t) — ; oC for ft 0 . 2 1 
|\0 for +40: | o(t) at = 1. 
2n-+] 
Hence b= — 0, a? a R.. am" — 0 (m0). (41 
47a“ 


(vii) Concentrated ring load over a circle 0 = 0,, —7 <b <7 and (i), 
This is a limiting case obtained from (38) by letting 


lim R(p)27a2(cos 6,—cos 6,) = Ro( p). (42 
6:05 
h-« 


This gives (40) with total load R,(t). The load per unit length of the circle , 


is Ro(t)/27a sin A,. 


The case of impulsive concentrated ring load follows from this, ané 
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results in (42) with load per unit length R,/27asin@,, where Ry has the 
same meaning as in (vi). 
Cast (viii) Uniform tractions and (i). If we take 0, = = in (iii) we obtain the 
se of uniform time dependent tractions applied to the cavity. Taking the 





,| limit as 6, > 7 in (37) or carrying out the integration in (30) we find 
} 
bm 0, ae 2R(p), am 0 (m,n + 0). (43) 
i Components of displacements and stress follow from (24) and (25): 
] F Ex+1 P 
s= — ’ R(: ax) S$ . —_—_— eif(a 1 v”) dé, 
Qra® . i(yé)?+E+% 
A 
: i(yéx)2-+-é2 - 
“a | R(E/ax) —\% ry SPT" pigic-1—w dé, 
2rx* , uys)"-+-o-+-t 
(44) 
os ' a y2)\( £o2)\2__ £4 
—. | R(é/aa) (1 = 2 MEx)— $0 —* ite 1—w ge 
Wl Aire” u(yé)?+E+2 
c U T Ty d ° 6d Q, O4 C4, 
ndent x r/a, Yy t/a, y = B/2a. 
rrving $s ; ve : : 
Using the convolution theorem and noticing that the inverse Fourier 
transform of R(€/ax) isaxR(aay), the first integral in (44) may be evaluated. 
{( To this end we note that the coefficient of R(é/ax) has two simple poles, 
both of which are located in the lower half of the é-plane. Hence a contour 
write | consisting of the real axis and a half circle in the lower half plane containing 
netion - both poles with centre at the origin may be chosen to determine the inverse 
Fourier transforms of the fraction in the integrand of wu. The convolution 
theorem then gives wu. Using this result in (2) and (3) we obtain o, and ag. 
u hap—1(4y?—1)-! rea-?(—sx+7)R*(y+1—2), 
{ 
r1 Ri aa(y+-1l—a 
2(4) 1)-tx-3 re|d(ysx)*+- sx i|R*(y 1-J]—2), (45) 
4) | 5) ~ We LR analy | -2x) | 
aye 3 wo |, (9,,2 »>y\2 ae lia DP  / | . 
(4) 1)—22 re |a(2y 1)(sa) — §2 Lil R \y +-1—z), 
cir 





Ria x7 )é is(y+1—x—n dn 
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valid for y > x—1. For y < x—1 all quantities are zero. If the function 


R(t) is discontinuous at t = 0 we define R(t) at t = 0 by 


0 for ¢t O— 
R(t) =(4R for t=0 
R for ¢=04. 


In return, when any of w, o,, og is discontinuous at y-++-1—2x 0 Fourier’s 
theorem gives one-half of the values obtainable from (45) at y+-1—a = 9, 
Physically y = 
ahead of the wave front the stress and the displacement components must 
vanish, as the analysis predicts, while at the points behind the wave front 
they are given by (45). On substituting the values of x and y from (44) we 
express the wave front by t/a = r—a. Here 1/a is the velocity of the 
dilatational waves. In this particular case, we only have the dilatational 
waves. This of course is also indicated by the type of boundary loading 
on the surface r = a. 

(ix) Uniform pulsating load and (i). 
be expressed by 


R(6, 4, t) R,e-. (46 
Hence R(0,¢,p) = 27R, 8(w—p), 
oe == 6, a} = 47R, d(w—p), a*=0 (m,n + 0). (47 
In obtaining (47) we used the formal relation 
© 
2775( p) ( ett dt 
Substituting the value of R given by (47) into (44) we obtain 
4uu/aR, = «2-*(Qar+i)[i(yQ)?—Q—7i]- exp[iQ(a«—1—y)], 
o,/Ry = x {i(yQa)?—Qe—i][i(yQ)?-Q—i]exp[iQ(x—1—y)], (48 
20/Ro c—{i(2y?—1)(Qax)?+-Ox+-t][i(yQ)?—Q- i|exp[iQ(a—1 —y)}. 
v= W = 7,9 = Tra = TH = O, Q = aaw, y>ar—l. 
(x) Blast loading and (i). In this case we have R -R, 8(t). This is 
a special case of (viii) in which R(p) = —Ry. The substitution of 


R(aan) 


- —R,d(axn) = 


— Ry 8()/ax 





x—1 represents the wave front. Naturally at those points | 


This type of normal tractions may 





ELAS' 
in (45 


2x 


valid 
the w 
substi 
we ha 


As inc 
u(1. y 

Fin 
spher 
can b 
volvii 


for th 


5. Ce 

Th 
comp 
The ¢ 
result 
for v: 
and f 
of th 
of the 
be no 











| ELASTO-DYNAMIC PROBLEM CONCERNING SPHERICAL CAVITY 267 


ction | in (45) gives the displacement and stress components, 
9 ] 
Qu 1 ; 1) u/Ry = x-*{(1—Tx)sin{ (2? —T*)*(y+1—2)]+ 
r(2' —T?)teos{(2P —T?)'(y+ 1—z)jeTwt1—, 
| 
9 1/2 i . 
f 1}? aao,/2Ro te 1}Px 1§(y+1—a)- 
riers t-a-86[ (1 —T)a?+-Ta—1]sin[(2P—V?)'(y+-1—2)|+ (49) 
1) 
yoint 1 (21°'—T"*)# (x? x)cos| (27° — "yay |] — x) |} Tw+1-z), 
must [2 \ , 2 1 ‘ 1s 
(=—1)* aaog/ Ry | 1}° (1- P)a—18(y+- 1 —ar) + 
front} = \I | 
4) y " cece . WP on 19\4 
; 2(1—T)?a2?— Tv 1]sin{ (2P—T?)4(y+1—a)+ 
of the 
tion (2° —T?)![2(1—T’)z?+-2]eos| (22° —T?)#(y+-1—2) Jje-Pot!-®, 
ding ' , ; ‘. - 2s Pe 7 —_— 
mas 4 O@> ( u ‘rv ‘rb "Od 0, I 27 ? 
valid for y > x—1. For the case of y < x—1 all quantities are zero. At 
the wave front y+-1—a = 0 one-half of the values obtained from (49) by 
(4g) substitution of y+-1—-2 0 should be understood. For example at ¢ 0 
we have ra a 0 
u(1, 0) s[p(A+2n)|-?*R, for t= 0 (50) 
47 p(A- 2) | *R, for ft @4.. 
\s indicated by Fig. 2 the value of u at x = 1, t = 0+ is the one to which 
l.y converges as ¥ >, 
Finally, it may be of interest to note that the problems dealing with 
spherically symmetric type of loading such as (i), (vi), (viii), (ix), and (x) 
n be treated more simply by solving the one-dimensional problem in- 
volving a single partial differential equation for wu alone. This was done 
lor the cases (viii) and (x) and the results were the same as in this paper. 
5. Computation and discussion 
48 The resulting formulae (49) in the case of blast type of loading were used to 
mpute the radial displacement wu and the stress components o, and 9. 
the computations were carried out on an electronic digital computer. The 
results are plotted on Figs. 2—4 versus the non-dimensional radial distance 2 
lor values of non-dimensional time y = }, 3, 1, 2, 3...., 11 as a parameter 
This is | “24 tor two distinct Poisson’s ratios v = } and }. Here y = t/aais the ratio 


ithe distance that the dilatational waves travel in time ¢, to the radius a 
the sphere. Some interesting results which follow from these figures may 


benoted. Att = 0a blast of impulse R, is applied to 2 = 1, the surface of 
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Fic. 2. Radial displacement produced by a blast. uy Ro/2u0(2/T 1)? 


Radial Stress 
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Radial Distance X 
Fic. 3. Radial stress produced by a blast. o,, = 2R,/aa(2/T 1)? 


the cavity. Immediately thereafter we find that a positive displacement 
u = Uy takes place. As time passes the displacement waves move into the 
media with sharp fronts ahead of which everything is zero. The variation 
of displacements in space and in time follows the interesting pattern shown 
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‘e notice, however, that after a while the tail of the wave front 


becomes negative, thus indicating negative displacements at some points 


behind but 

of the hoop 
lapse of tim 
behind it. . 
and then te 


Stress 
Oo 
O « 


oo 


4 
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becoming positive again later at the same points. The pattern 
stress og (Fig. 4) is similar to this except that after a sufficient 
e we find the maximum og not at the wave front but slightly 
\t the beginning oy is tensile, but later it becomes compressive 
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Fic. 4. Hoop stress produced by a blast 


the most unexpected result is that the radial stress in the 


od of the spherical cavity immediately after the application 


tensile. It starts with the zero value at x 1 and grows into a 


which is slightly behind the wave front. The mathematical 
he wave front y-+-1—a = 0 should be interpreted as explained 
the wave front we find infinite compressive stress. The study 
f rectangular pulse in the limiting case also indicates this fact 
more than the propagating blast pulse of finite impulse applied 
ice of the cay ity. 
ris indebted to Mr. J. C. Samuels for checking the analysis and 
out the computations. 
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THE INDENTATION OF A THICK SHEET OF ELASTIC 
WATERIAL BY A RIGID CYLINDER 


By DAPHNE G. PADFIELD and JEAN SIDA 


(Wool Industries Research Association) 
[Received 15 November 1956] 


SUMMARY 
The indentation of a uniform elastic sheet by a flat-ended cylinder has been com- 


ited for different ratios of cylinder radius to depth of sheet. The method used is 


direction application of the Hankel transform method described by Sneddon (1). 


1. Description of problem and method 
Tur indentation of the plane surface of a semi-infinite elastic medium by 
} arigid punch in the shape of a flat-ended cylinder has been fully discussed 
by Sneddon (1, p. 458). In this note we give an approximate solution of 
the corresponding problem for a thick elastic sheet. We envisage a uniform 
elastic sheet resting on a smooth rigid horizontal table, and wish to find the 
depth and the shape of the indentation caused when a heavy flat-ended 
circular cylinder is placed on top of the sheet (Fig. 1 (i)). The same analysis 
will also apply to the situation shown in Fig. 1 (ii) where the free sheet of 
material 1s acted upon by opposed loads on symmetrically placed rigid 
circular faces. The appropriate boundary conditions on the indented surface 
re: the depression constant within the circle of contact with the cylinder 
face, and the surface pressure zero outside this circle. 

From the theory given in (1), pp. 468 ff. it is easy to calculate the surface 
depression for any known pressure distribution, the numerical work being 
light when the function giving the pressure distribution has a known 
malytical Hankel transform. Three distributions of pressure p(r) over a 
‘ircular area (i.e. zero outside the circle) are particularly convenient; these 
are distributions proportional to 


i) l (0O<r<a) 
ii) a?—? (0<r<a) }, ef. (1), p. 528, 
iii) (a—r?2)-* (O0O<r<a) 


where a is the radius of the circle. For each of these pressure distributions 


x 


The Hankel transform p(r) is defined by 7(é) j rp(r)J,(ér) dr. The inverse transform 
0 


io @) 
ition (assuming convergence, etc.) is symmetrical, viz. p(r) { €D(€)J,(€r) dg. In the 
0 


work all transforms are for n 0. 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 3 (1957)]. 
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we may readily work out the corresponding depression. Thus the depression 
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Fic. 1. (i) Indentation of thick sheet on table by heavy cylinder. 
(ii) Indentation of free sheet by opposed forces on rigid circular faces. 
where d is the depth of the sheet, and 7(€) takes the forms 
(i) adJ,(ag)/é, 
(ii) 4aJ,(ag)/€—2a2J,(ag)/22, 
(iii) sin(aé)/é, 
respectively. The symbols v, 2 in equation (1) denote Poisson’s ratio an 
Young’s modulus. 
By taking a linear combination of the three cases (i), (ii), and (il 
u(r) = Au,(r)+ Bu,(r)+Cu,(r), say, we can approximate to the require 
boundary conditions, viz. 


u(r) = constant = u, say (O<r<a), 


surface pressure = 0 (a <1). 


In the actual computations three values of a/d were considered, vi 
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'ESSIOI 7 


ald = 1, 2,4. These results were supplemented by special solutions for the 

extreme cases a/d = 0,00 (see section 2). Foreach value ofa/d the depression 
(1)! u(r) was evaluated numerically from equation (1) for r/a = 0-0, 0-2, 0-4, 
(:6, 08, 1-0 as well as for a number of values greater than unity. Constants 
A. B. C were then chosen to minimize 


> w(r){Au,(r) + Bu,(r)+ Cug(r)—1}?, 


the summation being over the values r = 0, 0-2a, 0-4a,..., a, and the 
weighting factors w(r) being proportional to the areas of the corresponding 
ones. Knowing A, B, and C the ratio of the depression to the total load 
| was calculated, and the profile of the sheet for r > a was found. 


2. The extreme cases a/d = 0, © 

The necessary information about the case a/d = 0 is contained in 
Sneddon’s solution of the problem of indentation of a semi-infinite medium 
by a rigid cylinder. His result is that a total load 2#a/(1—v?) is required 
to produce a depression u - 

The case a/d = co can be examined by letting d > 0 in the expression 
for u(r) in equation (1). Thus, 


u(r 1—v2) —y 
r ! a Ep(é) (er) dé =“ a 


« 


0 


This shows that for 0 < r < a, p(r) = constant = (1—v?)-1Hu/d. 


3. Results 

Let ¢ denote the mean pressure applied by the circular base of the 
cylinder = total weight~za?. Then (u/d)/¢ = (1—v*)h/E where, from 
2, we find that h 0 for a/d = 0, and h = 1 for a/d = «©. The 





section 


results of the calculations described in section 1 show that 
h 0-716 fora/d = 1, 
h 0-844 for a/d = 2, 
h 0-903 for a/d = 4. 


The form of Sneddon’s result also indicates that h ~ tra/d ata/d = 0, i.e. 
it shows that the slope of the h v. (a/d) curve at the origin is }z. This curve 
is consequentiy fairly well defined for all a/d and is shown in Fig. 2. The 
profiles of the depression for a/d = 1, 2, 4 are shown in Fig. 3. 

The residues in the least squares calculations give an underestimate of 
ed, viz. the magnitude of the departure of the calculated depression from u = 1 
mae T 
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3. 3. Profile of depression for a/d = 1, 2, 
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The actual values of the coefficients A, B, C may also be of interest, since 
they determine the pressure distribution over the region of contact. The 


values were (in units of }#/(1—v?)): 





) —— ——_—_— 
79 7523 3°4921 
54 o'8211 0°3446 
57 1*5840 1°9879 
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EXPERIMENTS ON THE SLOW SWIRLING FLOW OF 4 | 


VISCOUS LIQUID THROUGH A TUBE 
By A. M. BINNIE (Engineering Laboratory, Cambridge) 


[Received 13 September 1956] 


SUMMARY 

To investigate the nature of the flow pattern when swirling water passes slowly 
through a long, straight, and transparent tube, a short length at the middle of the 
tube was caused to revolve, and through numerous holes in this portion the water 
was admitted under gravity. Colour injection was employed, and just clear of the 
rotating length three alternative régimes were seen: (I) downstream over the whole 
cross-section; (II) upstream near the axis and downstream near the wall; (III) 
downstream near the axis and the wall and upstream in the intermediate region. 
For four different arrangements of blockage in the two fixed lengths of tube, the 
results are plotted on diagrams having for their axes two dimensionless quantities 
that depend respectively on the tangential velocity at admission and on the mean 
axial velocity through the tube. When one of the fixed lengths was blocked close 
to the revolving part and colour was injected through the holes nearest to the 
blockage, the result under some conditions was three streams, symmetrically disposed 
120° apart, that passed through the other fixed length. 


1. Introduction 

WHEN swirling water passes through a convergent nozzle the axial com- 
ponent of velocity in part of the nozzle is sometimes reversed, and the water 
in that region returns into the reservoir from which it came. This effect, 
which was described by Binnie and Teare (1), clearly requires further 


_— 


— 


. . . . . *,° . ? 
investigation. It was obtained under complicated conditions. An air core 
existed, round which a boundary layer formed, having solid body rotation | 


and an axial velocity higher than the average over the whole cross-section. 
Moreover, a reservoir was essential in which the water from the external 
supply pipes could be steadied; consequently the velocity distribution at 
the nozzle entrance was unknown. 

To simplify the problem as far as possible, the nozzle was replaced by 
long transparent tube of uniform diameter. A short length of the tube was 
caused to revolve about its axis, and water was admitted through holes it 
this portion, thus the swirl (tangential velocity) at entrance was known, 
To conform with the theoretical case that Miss G. Vaisey is examining by 


s 
means of relaxation methods, the apparatus was made symmetrical, the 


perforated drum being at the middle of the tube. The velocities were kept 
very low, hence turbulence was avoided and an air core never formed. 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 3 (1957)] 
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In an arrangement of this kind, at least three flow régimes can occur in 
the fixed tube near the drum: 

(I) With the drum stationary, the water moves downstream at every 
point in the cross-section. 

(II) At the other extreme, with rotation of the drum but no discharge, 
a circulation is set up dividing the cross-section into two parts. Over the 
central zone the axial velocity is upstream, and downstream motion is 
confined to the annular region extending to the wall. It should be noted 
that the terms upstream and downstream applied to neighbouring currents 
do not imply a large shearing force between them, for both possess swirling 
motion in the same angular direction. 

(III) Under certain conditions a treble division occurs; the water near 
the walls and in the region of the axis moves downstream while upstream 
flow takes place in the intermediate zone. This effect, which was described 
qualitatively by Nuttall (2), appears to resemble reverse flow in a nozzle. 

The ranges of régimes I and II are not as small as has been indicated above. 
To find the limits of their existence and the conditions under which régime 
III occurs, four sets of experiments were carried out in turn with different 
arrangements of blockage in the fixed tubing: 

(A) Both the fixed lengths of tube were left unobstructed. 

(B) The outlet of the left-hand tube was blocked with a rubber bung. 

C) A sliding bung with a central hole was placed in both tubes at a 
distance of four diameters from the drum. 

(D 


The left-hand tube was blocked by a sliding bung fixed as near the 
drum as possible, and in the right-hand tube a sliding bung with a 
central hole was placed at a distance of four diameters from the 
drum. 

The numerical results are presented in the form of values of two non- 
dimensional parameters, which are based respectively on the mean axial 
velocity in the tube and on the swirl at admission. These parameters are 
defined by aW aV 

and =—, 


V Vv 


R - 


where a is the pipe radius, v is the kinematic viscosity, W is the mean axial 
velocity of the water in the tube, and V is the swirl velocity imposed on it 
at admission. The first parameter corresponds to the Reynolds number of 


flow without swirl. 


2. Description of apparatus 
Perspex tubing of internal diameter 2 in. was employed, its total length 
being 66 in., the central 4 in. of which formed the revolving drum. The 
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tubing was supported horizontally, as shown in Fig. 1, in an open tank 10 in, 
wide and 18 in, deep, which was divided by two diaphragms into a central 
compartment 50 in. long and two end compartments each 23 in. long. The 
drum was held in place by three wheels at each end, and water tightness 
between it and the adjoining fixed lengths was maintained by Gaco rings 
lubricated by means of screw-down grease cups (not shown). It was rotated 
by an electric motor and a variable-speed gear box placed outside one end 
of the tank, the drive being transmitted through a long horizontal shaft toa 
small gear-wheel which engaged with a toothed ring on the drum. Water 
was led under gravity from a constant-level tank to a submerged inlet in 
the central compartment. After passing into the drum and dividing into 
equal streams in the two fixed lengths of tube, it reached the end compart- 
ments, which were connected by external piping to a small tank and weir, 
The flow through the apparatus was controlled by altering the heights of 
the constant-level tank and of the weir, and it was easily measured with 
a graduated flask held under the weir for a known time. The tests were 
usually carried out with the tube in the central compartment submerged 
a couple of inches, and even at the highest discharges the outlets of the fixed 
tubes were well covered. The perforations in the drum consisted of 11 
belts, each of 60 holes 0-080 in. in diameter, symmetrically disposed 0-25 in. 
apart. The total cross-section of the holes was about 5 per cent greater than 
the cross-section of the tube. The wall thickness of the drum was 0-25 in., 
so that it is reasonable to take the swirl velocity imposed on the water at 
admission as equal to the tangential velocity of the inner surface of the 
drum. The revolutions were timed with a stop-watch, and the maximum 
speed that could be obtained without undesirable vibration was about 
40 rev/min. 

In order that the flow pattern inside the tube might be examined, three 
tappings designated R,, R,, R, were drilled along the top of the right-hand 
fixed tube at distances 2-65, 5-65, 8-65 in. from the central plane of the 
drum, and likewise three tappings L,, L,, L,, were added symmetrically 
to the left-hand tube. At the tappings, coloured liquid was slowly injected 
through a probe of outside diameter 0-043 in. mounted on a vertical catheto- 
meter. The tip of the probe could be traversed across the diameter of the 
tube but, to avoid unnecessary disturbance, only the upper radius was 
usually explored. The coloured liquid was a solution in water of the sodium 
salt of eosin, made up to density 1-003 and then diluted with methy] alcohol 
to density unity. 

The apparatus thus differed from that used by Talbot (3), who employed 
a long unperforated tube in which swirl was induced by the rotation of 
part of the tube situated some distance from the entry. 
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Side elevation of apparatus 


Fia. 
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3. Experiments on symmetrical flow with no obstructions in the 
fixed tubes (arrangement A) 


Preliminary tests first with no holes in the drum and then with a single 
row of 30 holes 0-020 in. in diameter demonstrated the existence of régime II 
and the necessity for a great number of holes in order to avoid turbulence 
at entry. Accordingly the 660 holes mentioned in section 2 were drilled; 
for reference the eleven rows of holes will be numbered from —5 through 
0 to +5 reading from the left. To check the symmetry of the flow a bung 
was inserted at each outlet of the fixed tubes in turn, and with the outlet 
weir unaltered the discharge through the two halves of the apparatus was 
found to be the same. Colour was then inserted into the middle compart- 
ment outside but very close to the central row of holes, and after admission 
it was seen to divide symmetrically into two streams. The experiments of 
Binnie and Teare (1) with a double-ended swirl chamber operated at high 
swirls showed large axial velocities across the central plane, and a similar 
effect in the present apparatus was thought possible. But it was only on 
rare occasions that colour injected at R, found its way into the left-hand 
tube. In these few instances colour under the influence of régime II was 
carried up near the axis to the central plane where it remained visible for 
a matter of minutes, and occasionally a small proportion of the colour 
drifted from this stagnant zone across the central plane. 

The usual method of test was to set the discharge and then to increase 
the speed of rotation in steps. At each step a measurement of the discharge 
was made which, owing to centrifugal action, decreased slightly with rising 
speed, and at R, colour was injected at r/a = 0-9, 0-8,..., 0, where r is the 
radius at which the tip of the probe was placed. The kind of régime was 
noted and the result plotted in Fig. 2. The two sets of observations at 

values of R about 140 and 260 were made when the outlet of each tube was 
obstructed by a bung with a }-in. hole; the results are seen to be in agreement 








) at ria 
ties; b 


. . . . ‘a } 
with the tests under normal conditions, and this was found to be true also 


of a set at about R = 76. At V = 0 and the largest values of R, upstream 
flow was found very close to the wall. This was part of a local eddy induced 
by the high velocity of the water at admission. It was thought reasonable 
to plot this observation as régime I, which was found at all lower values of 
R with the drum stationary. 

To explain the nature of régime IT an account will be given of what was 
seen at V = 460 (8 rev/min), R = 104 (W = 11-3 in./min). Under these 
conditions the coloured streamline from a fixed position of the probe did 
not constantly follow exactly the same path, but its wanderings did not 
obscure the general flow pattern. With injection at R,, the colour inserted 
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tra = 0-9 was carried downstream with vigorous swirl and axial veloci- 
ties: both components greatly diminished beyond R, but the swirl was still 
pe reeptible as the colour entered the end compartment. It should be noted 
that inevitably the colour left the probe w ith a small radial velocity, there- 

fore the effective radii of injection were somewhat less than the values given 
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Fic. 2. Régime diagram for arrangement A 
Régime I 
e 9 II 
here. At r/a = 0-8 and 0-7 the initial pitch of the helix was smaller than 
before, but it rapidly increased downstream as the swirl died away. At 
ra = 0-6 the initial axial velocity was zero; some of the colour was whirled 
way downstream and the rest moved upstream shedding filaments which, 
iiter moving to a larger radius, were caught in the fast downstream flow. 
Some succeeded in moving as far as row +-5 before being turned back. At 
smaller values of v/a all the colour initially moved upstream, the shedding 
process in the fixed tube becoming less marked. At r/a = 0-3 all the colour 
entered the drum, in which its swirl and upstream velocities were very small 
intilit drifted to a larger radius and was suddenly caught in the main stream. 
At values of r/a < 0-3 the colour reached the central plane but with extreme 
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slowness in the last part of its travel; from there a trace escaped into the | 
left-hand tube. At r/a = 0 the colour on the central plane concentrated 
into a disk having a radius nearly half that of the tube. From the periphery 
traces were carried away downstream by the main flow, and a small propor: | 
tion was removed into the left-hand tube in a similar manner. 

These observations were confirmed by a traverse at R,. The position of 
zero axial velocity was at about r/a = 0-45, showing that the area of the 


core of upstream flow was diminished by a half in the length between R, 
and R,. External injections were also made in close proximity to rows 0, 
2, and +5. From 0 the colour divided symmetrically and passed down- 
stream at a little distance from the wall; but when the drum was clear and 
a minute had elapsed after the end of the injection, traces returned simul. | 
taneously into the drum from both fixed tubes at a radius of about } in, 
From row +2 most went swirling away past R, but a little, after reaching 
a point between Ff, and R,, returned slowly at a smaller radius into the 
drum. From +5 the same reverse trace was seen, but the main flow was 
downstream so close to the wall that beyond R, most of its swirl was lost. 

Almost the whole of Fig. 2 is occupied by régime II. Along the horizontal 
axis of the diagram régime I is bound to exist, but a very small amount 
of rotation sufficed to remove it, consequently it is confined to a narrow 
band at the bottom. Régime III was not found within the ranges of V and 
R that were explored; and as this type of flow was the principal object of 
the investigation, an attempt, described in the next section, was made to 
produce it by destroying the symmetry of the apparatus. 


4. Experiments on asymmetrical flow produced by blocking the 

outlet of the left-hand tube (arrangement B) 

In this series of tests the outlet of the left-hand tube was closed by a 
rubber bung. The alteration had the desired result, for over a certain range 
of R shown in Fig. 3, régime I changed to III when the speed of rotation was 
sufficiently increased. With the usual procedure of injection at R,, an 
almost constant discharge, and a start with zero rotation, the decay of 
régime I as the speed was raised was indicated by a fall in the axial velocity 
at about half radius, and the velocity at this position was eventually re- 
versed. As R was increased, a higher rotational speed was required to effect 
the transformation. Fig. 3 differs greatly from Fig. 2. The observations 
along the axes were the same as before, but régime II appears nowhere else 
on the diagram, for it was eliminated by a slow rotation and replaced by I 
at low axial velocities and by III at high. 

More complete observations were made of régime III at V = 1,79 
(28 rev/min), R = 68 (W = 6-6 in./min), under which conditions the 
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streamlines were steadier than in the case detailed in section 3. With R, 
in use, the colour injected at r/a = 0-9 and 0-8 was carried downstream 
out of the central compartment, losing its swirl before it had travelled far. 
Atr/a = 0-7 it circled upstream } in. and was then carried away at larger 
radius by the main downstream flow. This effect was more pronounced at 
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radii down to r/a = 0-4, the colour getting just inside the drum before being 
swept downstream. For smaller radii the initial axial velocity was down- 
stream, though at r/a = 0-3 the travel in that direction was slight, being 
succeeded by a very slow movement just into the drum at a larger radius. 
\t ra = 0-2 much the same was observed, but the colour drifted into a 
stagnant zone in the form of a conical frustum with its smaller end about 
, i. downstream from R, and its larger end in the drum between row +-5 
nd the nearby joint. At r/a = 0 the cone was completed with extreme 
slowness, its vertex lying on the axis about 1 in. from R,. Thus, in this 
example, with the rotational speed not far above the minimum at which 











284 A. M. BINNIE 





régime IIT could exist, the downstream velocity in the central part of the | at R 
) ] 


cross-section was feeble, but even when the speed was raised to 40 rev/min, | 


the length of downstream travel at the axis was not noticeably increased, 
With injection at R,, the motion at values of r/a down to 0-5 was down- 
stream with a slight swirl that was soon lost; at r/a = 0-5 some of the colour. 
after going an inch, returned at a slightly smaller radius and reached the | 


nearer joint, then being at a larger radius. At r/a = 0-4 and less, the motion } 


was upstream just into the drum at a radius greater than that of injection: 
evidently the flow was passing outside the cone mentioned above. 

The most convenient way of exploring the blocked half of the tubing 
was to continue to inject into the right-hand tube, the bung being trans. 
ferred to the outlet of that tube, but the observations will be recorded as 


if they had been made in the lef:-hand portion. At L, and r/a = 0-9, the 


colour vigorously swirled away from the drum, and after reaching J, it 
returned at a smaller radius and passed through the drum, losing most of 
its swirl as it did so. At r/a = 0-7 the flow was directly into the drum as 
far as row —5 where it was caught in the stream nearer the wall and turned 
back into the blocked tube; eventually a trace returned into the drum at 
a smaller radius. At r/a = 0-5—0-3 all the colour entered the drum; some 
got through to the right-hand tube at a larger radius with a velocity that 
was low compared with the brisk motion of the rest which was returned 
towards L,. At r/a = 0-2 and less, the whole of the colour passed through 
the drum and emerged in two streams, one moving slowly at a large radius 
and the other with almost zero swirl and axial velocities close to the axis. 
This striking effect was particularly well marked at r/a = 0. The colour 
all reached the central plane where most moved quite abruptly to a larger 
radius and entered the right-hand tube. After it had gone, the remainder 
could be seen near the axis; it advanced with extreme slowness and finally 
came to rest about 1 in. downstream from R,, which was the position of 
the apex of the cone described in the previous paragraph. All the examples 
of régime III that were closely examined showed this division of the stream 
into two very different parts. At Z, it was only extremely close to the wall 
that any flow existed away from the drum; the stream swirled towards the 
bung for 3} in. and then crossed the section at a much smaller radius. At 
r/a = 0-9-0-4 the colour moved initially without swirl towards both the 
axis and the drum, but later in its travel it picked up swirl and most of it 
entered the drum. At still smaller radii the colour moved towards the drum 
with swirl that was evidently derived from the flow near the wall. An 
attempt to inject colour externally was unsuccessful because the inward 
velocity was so slow that the colour diffused before it entered the drum. 
To discover to what extent the position of the blockage altered the régime 
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a sliding bung was placed in the left-hand tube, and two sets of 
experiments were made at constant V and R and varying positions of the 
blockage. In the first, at V = 1,420and R = 92, régime III was found when 
the blockage was 4, 84, and 16 in. from the nearer joint, and II when this 
distance was 0, 1, 2 in.; in a trial at 3 in. the axial velocity near the axis 
was impercept ible. In the second, at V = 2,360, R — 84, the limit between 


| the two régimes was at 4 in. Finally, the appearance of régime II when 


the blockage was near the drum was confirmed by fixing the bung as close 
to the joint as possible and making observations at constant R = 116 and 
varying rotational speed. Over the range 470—1,890 of V that was examined, 
oily régime II occurred. At the highest speed the motion within the drum 
it small radii was very slow, and colour injected at the axis, after approach- 
ing the bung, moved into a clearly marked zone in the form of a conical 
frustum with its smaller end close to the bung and its larger near row +5. 
From the larger end the colour was whirled away downstream. 


5, Experiments on symmetrical flow with a perforated bung in 

both the fixed tubes (arrangement C) 

Asliding bung with a central }-in. hole was now placed in both fixed tubes 
at a distance 8 in. from the nearer joint, thus the apparatus was converted 
intoa primitive kind of double-ended swirl chamber. The results are plotted 
in Fig. 4. Along the axes the observations at R, were the same as before, 
and at a very low value of R there was, with rising V, the same transforma- 
tion from régime I to II only, that is shown in Fig. 2. Under these con- 
ditions it seems that the bungs were too far away to have much influence 
on the flow at R,. At greater values of R, however, régime III was visible, 
the necessary rotational speed for its appearance becoming less as R was 
increased. At first, at R = 60 approximately, régime ITI was of the feeble 
kind described in section 4, in which the downstream travel over the central 
zone was small. But above R = 200, in all the examples of régime III that 
were observed, the central flow consisted of a vigorous and clearly defined 
stream that passed through the hole in the bung. 

At V = 1,990, R = 260, with external injection at row 0, the colour 
divided symmetrically and soon formed into a diffuse cylinder at half 
radius extending over the length of the drum; the cylinder faded as colour 
was lost downstream from its ends, but none moved inwards to supply the 
discharge along the axis. At rows +-2and +5 the colour was rapidly whirled 
away downstream at a large radius. With injection at R,, the colour at the 
largest radii moved downstream and was quickly diffused, but as r/a was 
decreased below 0-7, reversed flow made its appearance. At r/a = 0-5 all 
the colour started upstream; some of it was removed downstream at a larger 
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radius, but the rest entered the drum in a rather diffused state and later 
formed itself into a cylinder of small radius extending from row -+-5 to R 
and persisting for several seconds. At r/a = 0-3 the initial axial] velocity 
was slight, nevertheless the cylinder formed again, extending quickly down. 
stream to the hole and slowly up into the drum in a more diffuse state, 4; 
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Fic. 4. Régime diagram for arrangement C 
Régime I 
° 99 I] 
II 
smaller radii the motion was entirely downstream in the form of a sharply 
defined helix of small pitch that disappeared through the hole. Occasionally 
the filament bounced back out of the hole at a slightly greater radius, 
indicating that reverse flow existed there also. 


6. Experiments on asymmetrical flow produced by blocking the 
left-hand tube close to the drum and inserting a perforated bung 
in the right-hand tube (arrangement D) 

Finally a single-ended swirl chamber was simulated by placing in the 
left-hand tube a sliding unperforated bung as close as possible to the drum 


and in the right-hand tube a sliding bung with a }-in. hole at a distance of | 
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gin. from the drum. The results are plotted in Fig. 5, which shows that 


wain régime IIT was readily obtained. As in Fig. 4, for a rather small value 
fR régimes I, II, and III appeared in turn as the speed was increased, the 
nset of the change from II to III being again indicated by the decay and 
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eventually the reversal of the upstream velocity at the axis. At higher 
values of R the intermediate régime II was not observed. 

The details follow of what was seen of régime III at V 1,080 
min), R 134 (i 12-4in./min). With injection at R, the colour 
inserted at r/a = 0-9 went downstream at a high axial velocity and diffused 
near R,. At r/a 
of the colour came back at a smaller radius and entered the drum. The 


l6 rey 
0-8 after moving downstream a couple of inches, some 


same was seen at r/a 0-7 except that a portion of the band that reached 


the drum was carried briskly downstream at a larger radius. Between 
ra = 0-6 and 0:4 all moved towards the drum and some remained stagnant 


inside it; the remainder was removed downstream at a larger radius when 
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r/a = 0-6 and at a smaller radius when r/a = 0-4. The motion at r/a = 03 
was in both directions, part going very slowly to the hole and the rest far 
into the drum at an increasing radius. For smaller values of r/a the colow 
moved downstream with its axial velocity increasing as it approached the 
hole. 

With external injection at row +5 the colour swirled rapidly downstream 
close to the wall and diffused beyond R,. At row 0 the downstream axial 
velocity was at first so slight that the helix of colour that formed was barely 
distinguishable from a cylinder. At row —5 the colour moved into contact 
with the face of the bung where it could be seen in-a diffuse state. Then, 
after a pause of several seconds, it gathered itself into three sharply defined 
streams 120° apart at a radius of about 3 in., which passed with considerable 
angular velocity through the drum and far down the right-hand tube. A 
horizontal photograph, taken with a vertical camera and a submerged 
mirror inclined at 45°, is reproduced in Fig. 6, which shows the three streams 
with their heads almost in line when their development was well under way. 
The effect was found to be independent within wide limits of the duration 





and intensity of the injection, but one of the streams was faint when the 
duration was appreciably less than the time of one revolution of the drum; 
no alteration in the effect was noticed when simultaneous injection was made 
from a second probe fixed over the same row and distant about seven holes 
(42°) from the first probe. Similar results were observed when injections 
were made at rows —4, —3, and —2. The number of streams seen in many | 
repetitions of the experiment was three, but four occurred on a few occa- 
sions when, apparently, conditions had not become steady after alterations 
in speed and discharge. 

To elucidate this matter further, the discharge was left unchanged and 
observations were made at various speeds with injection at row —5. At 
8 rev/min (régime I) the colour again collected up against the bung, and 
after a pause a helix of increasing diameter emerged giving rise to a cone 
of colour that passed downstream, though traces long remained visible in 
the drum. Thus the single point of injection gave rise to a single stream. 
The same occurred at 12 rev/min where régime III had just set in. At 
24 rev/min the three streams were not quite as sharply defined as at 16 rev 





min; at 32 rev/min they were definitely more diffuse, and whether they 
numbered three or four could not be decided. When this set of experiments 
was repeated with R increased to 216, it was found at the lower speeds that 
the cone tended to concentrate into streams and, as it passed through the 


a 


drum, it threw out streams ahead. However, at 32 rev/min three sharply 
defined streams emerged right from the start, as in Fig. 6. 
In conclusion, the perforated bung was removed, and thus the kind of | 
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blockage was restored that is considered at the end of section 4. With the 
game discharge and swirls that are mentioned there, tests were made with 
external injection at row —5. Under these conditions régime II had been 
ven at R,, and now inside the drum separate streams of colour again sprang 
from the neighbourhood of the bung. But they were much more irregular 
both in spacing and in length than when the perforated bung was in place, 
nd their number varied from two to four. 


7, Further discussion of results 

Figs. 2-5 differ widely, except along the axes, although the boundary 
mditions at entry remained unaltered; thus the arrangement of the 
blockage was of paramount importance. The chief object of the experiments 
was achieved in that they showed that régime IIT could be produced, but 
the velocities that were reached were small compared with those that occur 
nanozzle with an air core. It appears that viscous swirling motion in a 
tube is so complicated that little can be done to explain it in general terms, 
ut the existence of régime II at zero or small discharge can be accounted 


fry 
1 





by the centrifugal action inside the drum which caused the pressure near 
the wall to exceed that at the axis. Thus the water near the wall was forced 
terally into the fixed tubing, and at the axis a current in the reverse 
lirection was set up. Régime III is harder to explain. In the vigorous form 
sen in arrangements C and PD the existence of the holes in the bungs 
stimulated a strong downstream flow close to the axis. Now it is well known 
that on the back wall of a single-ended swirl chamber there is inward flow 
cause the liquid is retarded near the wall, its swirl being reduced so that 
| itcannot maintain its position in the centrifugal field of force. Hence over 
the faces of the perforated bungs an inward stream was present, and it seems 
that in the narrowing region between this and the stream close to the axis 
some of the water was forced to move upstream. The origin of the central 

lownstream current varied with the arrangement of the blockage. In B 

the source was the blocked left-hand tube in which régime II existed, and 

Cand D the current seemed to arise from the inner surface of the reverse 

mnular stream. In A, where régime IIT was not seen, there were no per- 

lorated bungs to enhance the downstream motion on the axis nor a blocked 

| tube in which it might develop; and although it was looked for, inward 
\ radial flow from row 0 remained undetected. 

The effect shown in Fig. 6 does not seem to have been previously noticed. 
For the reason set out above, the colour inserted outside row —5 moved at 
irst to the face of the neighbouring bung, but its emergence in separate 
‘teams was unexpected. The passage of these streams past A, accounts 
i part for the unsteadiness of the colour bands injected there. 

“ 
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\ JET DEFLECTED FROM THE LOWER SURFACE OF 
\N AEROFOIL 
By H. J. DAVIES and A. J. ROSS (University of Southampton) 
Received 26 July 1956] 


SUMMARY 


sional, steady motion of an infinite, perfect, incompressible fluid 
1 thir mmetrical aerofoil at zero incidence to the mainstream is considered 

t ies from a general point on the lower surface of the aerofoil at an angle 

hord. The motion in both the mainstream and jet is taken to be irrotational. 
letermination of the mainstream flow about the aerofoil and jet, the aerofoil 
presented by a single plane vortex sheet and the jet by a single curvilinear vortex 
TI trength of the latter is derived by considering the flow within the jet. 

An iterat method of solution is suggested for the evaluation of the true jet shape 
re mn the aerofoil; the first step requires the assumption of an approximate 


\s an example, the lift, pressure distribution, moment coefficient, and 
| | 





f pi ive are calculated for a jet exit position close to the trailing-edge with 


if ejection and momentum coefficient. 


1. Introduction 


We shall consider the two-dimensional, irrotational, steady motion of an 
finite, perfect, incompressible fluid past a thin, symmetrical aerofoil at 
ero incidence to the mainstream when a jet issues from a general point 

| on the lower surface of the aerofoil at an angle to the chord. The fluid 
| motion in the jet is assumed to be irrotational and bounded by curvilinear 
rtex sheets across which the pressure is continuous, and the velocity 
uscontinuous. At infinity the jet will be parallel to the mainstream, as 
lollows at once from momentum considerations. The condition that the 
uid is perfect implies no mixing of the jet with the mainstream which, in 
turn, implies constant total pressure heads H and J in the mainstream and 
t respectively. Under these circumstances, the total lift experienced by 
 aerotoil, i.e. the pressure lift plus the component of the jet reaction 
ibe given by pl’ ¥ (+) where p is the density of the fluid, U the 
unstream velocity, and 5 (['+y) the total circulation around a contour 
| at infinity (1, 2). Hence, it is required to find the sum of the vorticity on 
the aerofoil $ and on the jet boundaries > y. 


nm 


. The mathematical model 

Suppose a jet of small cross-section B, B, issues at an angle 7, from the 
wer surface of a thin symmetrical aerofoil which is at zero incidence to 
the mainstream U as shown in Fig. 1. The aerofoil is replaced by a vortex 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 3 (1957)] 
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sheet of strength I’ per unit length measured clockwise, and the jet boun-| 
daries B,C,, B,C, by vortex sheets of strengths y, and y, per unit ar 











length, also measured clockwise. 
| 
y 
(a) 
U { 
(b) 





Fic. 1. (a) The physical plane, (b) jet section 


Consider the one-dimensional flow across any section d of the jet across 
which the mean velocity is V and the mean pressure p. Subscripts 0 and«| 
will be used to denote conditions of the jet at the exit and infinity. As B, 
and B, are stagnation points in the mainstream the pressure at the exit of 
the jet will be equal to the total head H of the mainstream. Hence 


4p(V2—V2) = H—-p>0, ie. V >}. 


Thus, in the absence of mixing, the jet boundaries converge from the exit 

to infinity. The square of the ratio of the jet cross-sections d, and d,,, from 

continuity considerations, is given by d?/d?, = V2,/V2, which, by applica 

tion of Bernoulli’s energy equation, may be written as d2/d2, = 14+U?/I 

Thus, it is seen that if the jet velocity V, is very much greater than the 

mainstream velocity U, the jet may be regarded as having a constant cros*| 
section. Finally, following Taylor (3) and Spence (4), the limit of the jet! 
width is taken in such a way that the mass flow tends to zero but the? 
momentum remains finite and constant. The jet is then represented by: 

single vortex sheet BC of strength y per unit arc length, across which thi 

pressure as well as the velocity is discontinuous. 


3. The vorticity on the jet boundary 

A section of the jet is shown in Fig. 1 (b), q and q’ being the velocity| 
outside and inside the jet respectively, p the pressure, and the subscript: 
1 and 2 referring to the boundaries B,C, and B,C, respectively. 

The strengths per unit arc length of the vortex sheets B,C, and By“ 


regarded as boundaries of the flow in the mainstream will be y, = - Iv} 


Yo = 4s. If corresponding increments of are lengths of B,C,, B,C, a0"! 
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BC are taken to be equal the vortex density on BC will be given by 
y=¥,+Y2 = %—G- As the pressure across B,C, and B,C, is continuous 
and the total heads of the jet and mainstream are constant an application 
of Bernoulli’s energy equation across the section shown in Fig. 1 (6) gives 


G—H = %2"°—G;" = 2(p,—Pe2)/p 


so that 1o—h f . 2, where V = 3(q;+4)) 
pP 
d grad ; , 
va since d is assumed small, 
ph 
dV? 
VR 


by consideration of the acceleration normal to the jet, R being the radius 
f curvature of BC. Hence 7 


Vd 
y = we (3.1) 
UR 
where | 4(q, +d») is taken to be constant and equal to the mainstream 


velocity U’. It has already been stated that the model considered is one in 
which the momentum of the jet remains finite and constant. So that 
pl*d = py Vido, where po, the density of the jet fluid, is taken to be equal 
0 p. By defining the non-dimensional jet momentum coefficient as 
Cy = po Vod,/4pU 7c, where c AD the chord length, (3.1) can then be 
written as 


LUC, ¢/R. (3.2) 
4. The vorticity on the aerofoil surface 

Consider rectangular axes in which Ox is parallel to and in the direction 
of the mainstream and Oy is normal to the aerofoil chord AD, with origin 
it the leading edge A (see Fig. 1). One of the boundary conditions of the 
flow in the mainstream is that the velocity normal to the aerofoil surface 


is zero. From consideration of the model described in section 2 this con- 


dition gives 
l fT(l)dl yv(a—X) ds 
i | Wea? @O<e<d, (41) 
a7) «—l 2r | (a—X)?+Y? 
where / is the distance of a general point on the chord from the leading 


edge and X,Y the coordinates of a general point on BC corresponding to 

inare length s measured from B. The second integral in this equation may 
interpreted as the downwash on the aerofoil due to the vortex sheets 
nding the jet and will be denoted by —(1/27)w(2). 
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Consider the integral 


x(c—ax)! dx 
(A—a)(a—l) 


0 
which, by substituting x = csin?6, becomes 


i sin?@ cos?4 dé oe 
- See eae erat A= (4,2 
J {(A/c)—sin?6}}sin?6— (l/c)} 
0 
= * 23(c—2)} 
Thus the operator = dx 
(A—2x) 
0 
is seen to be a null transform of order one to the kernel | /(x , in the 
interval 0 < # < c. By the application of this transform to the equation 
yy 
) dl = w(x) 
| xl 
0 
it is seen that 
wee -  gi(c—2)} * wlx)at(e—ax 
| rm | YY  dadl+- R(A) = "de, (43 
; J (A—2x)(x—l) : A—2 
0 0 0 


where F(A) is the residual, arising from the singularity in the integran 


when the order of integration is reversed. Thus 


A € A € 
R(A) = lim} - rl) ih Ans ee 
e—0 P x (A r)\(xv l) 
l \ € I A-e« 
A+e 
» &i(c—2) PO and 
} Ne J 2 
r A—e€ l A—e€ 
which, on substituting / = e€+A, x = en+A and proceeding to the limit 
gives R(A) m72T(A)AL(e—A) (4.4 


(see reference 5). Thus the inversion of the integral equation (4.1) is, from 


(4.2) to (4.4), given by 


c 


- ] " F OL )ezE(C xr 
ra) = — z | Pl) dl— | Se te 
7"A (Cc A) A - 


where 


| 
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/ The constant | ['(/) dl is determined from the condition that B is a 
0 

stagnation point of the mainstream. The condition that the velocity 

| normal to AD is zero has been satisfied above. Hence, to satisfy the 


dition that B is a stagnation point, it is sufficient to write down the 


; , 
| condition that the 2-component of the velocity at that point, taken to be 
"i / h iS zero 1.e€ 
D(b l ' y) 
P el ds 0, (4.6) 
2 or (b—X)?+ Y? 
Be 
here P denotes the Cauchy principal value of the integral. This gives 
» avalue for [(b) which may be substituted in equation (4.5), writing « = b, 
rder to obtain the value of | [(/) dl. 
t is of interest to note that the significant difference between the above 
lel and that adopted for treatment by thin-aerofoil theory (4) is the 
position of the jet exit, as this changes the boundary condition required 
r the determination of the constant arising in the solution of the integral 
equation for the vorticity on the aerofoil. This results in a different order 
f infinity of the vorticity on the aerofoil at the trailing edge, and an 
nite difference in the vorticity on the jet mean line there, although 
orticity is unchanged. However, it can be shown that the 
when the jet is taken at the trailing edge as in reference (4) is the 
s the jet position tends to the trailing edge in the solution contained 
5. The lift 
he total circulation within the mainstream is given by 
(rajd+ [ yds, 
e 
Ls tated in section 1, the lift experienced by the aerofoil will be 
L ol} | (1) dl 3 | yds . (5.1) 
| 0 LB 
by substituting for y from (3.2) and transforming the variable of integra- 
n {roi c length to tangent slope so that ds R dr, the total vorticity 
the jet boundaries may be written as 


| yds = 30,7, Ue. (5.2) 
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The method of calculating the circulation around the aerofoil has been | 0n th 
given in section 4. By substituting for y and 4s, given above, in equations | ejectic 
(4.5) and (4.6), together with A = 6 in (4.5), and eliminating I'(b) between | mean 








these two equations, the aerofoil circulation is then given by the equation | angle: 
c 0 equati 
r on , €,0e f Y dr . 
T(l) dl = mb*(c—b)!| 2U ——2 — yey + — 
L 2a (6—X)?+Y? equat! 
0 To 
P In « 
ey r3(Cc—2 region 
So | ool) —*) dx, (5.3) |, ° 
* b—x be det 
0 
the as 
0 rt | 
, ' r r— Bf siope 
where w(x) = —43C,Uc ——. dr. (5.4) 
ial (a—X)?+Y? true p 
i ever, i 
It is seen from equation (5.3) that another boundary condition is required | of the 
to determine the remaining unknown quantity, namely the equation of the | ment. 


limit line BC of the jet. The outstanding condition is that the velocity of | at the 
the mainstream normal to BC is zero. By taking X*, Y* to be the co- | 
ordinates of a general point and r* the slope at that point of the line Bé 


given by X = X(r), Y Y(7), this condition may be written in the forn 
C ry} In thi 
i skas, 21 ff ()Y* 
sin r*| U 4 is ' fe wre p 
| 27 | | (AP 67-2 * then 
0 } } 
: r r% ; | 

(2 J c) . (X * Es )? | () * } )? ‘ bee 
To Th 
eS rp Tr rx ) leadit 

COS T* | D(A 1) dl | 

In J (X*#—12?4+-Y*" | 

0 ) 


(<*—Z) 


+(4C,Uc)P | —— dr}. 5.5 
(3 B ) (X*- Xp(y* yp | \ 


. 


It does not seem possible to obtain X and Y as explicit functions of 7. 


However, if an approximation to the equation of BC can be obtained, then | wher 


equations (5.1) to (5.3) will give a first solution for lift whilst equation 


(5.5) offers a means of iteration. 


6. The limit line 
As an approximation to the equation of the limit line BC of the jet it | Thy, 
was decided to use the solution of the source type flow given by Woods (6). | 
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On the assumption that jets of equal momentum and at equal angles of 
ejection have the same influence on the mainstream, the equation of the 
mean line of Woods’s source type flow, having the same momentum and 
ngle of ejection as that of the jet, was taken as an approximation to the 


| equation of BC. The downwash w(x) on the aerofoil may now be evaluated 


yds (6 





numerically and the circulation around the aerofoil obtained from 
equation (5.3). 
In order to obtain a new position of the jet limit line the flow in the 
region of the assumed curve must be considered. Now tan7*, which may 
e derived from equation (5.5), will be the slope of the streamlines cutting 
the assumed curve and must, with each successive iteration, tend to the 
slope of the limit line itself. The velocity field of the flow indicates the 
ie position of the jet but the computation involved is too lengthy. How- 
er, it is possible to derive a measure of the displacement along the normal 
the assumed curve from the normal velocity. If is the normal displace- 
ment, s the are length, g,, the normal velocity, and q, the tangential velocity 


the point on the assumed curve then dn/ds = q,,/q, i.e. 


Cn ds 
J 4 


aie 


it] 
this way a new equation to the limit line may be obtained, since if X,, Y; 


e points on the assumed curve and X,, Y, points on the corrected curve, 


x, = Z, Y, = ¥,4 


nsinr and 


7. COST. 
7, Moment coefficient and pressure distribution 


[he pressure force contribution to the moment coefficient about the 


ding edge of the aerofoil will be given by 


p Pu )x dx 


| C7 be Gy je aa 


ul » | (1+U,) P(a)ax de. 











298 H. J. DAVIES AND A. J. ROSS JET D! 


The pressure coefficient is given by C,, = 1—gq?, where q has the values ¢, By we 
and q, given above. ; so that 
8. Numerical analysis and 


For the numerical work it is more convenient to use a dimensionless form | yhere 
of the equations; lengths are given as ratios of the aerofoil chord and veloci- 
ties as ratios of the mainstream velocity at infinity. The equations of th 
preceding paragraphs will remain unaltered in form if ¢ and U are replaced | the lin 
by unity. 

The evaluation of the principal value of integrals involved in this work 


is mainly standard work. However, there occur two integrals the evaluation 





The in 
of which calls for an explanation. Consider the integrals 
0 om . 
| : - —. dr (8.1 ; 
J (XxX X)?+ ( i ' 
‘ | Simila 
. go | 
and | ia = Ar (8.9 
(X*—X )?+-() Y)? 


occurring in equation (5.5). Let P*(X*, Y*,7*) be any point of the cw 
BC, P(X,Y,7) a neighbouring point of are length s and chord length } | 9. Nu 


from P* and 6 the angle between P*P and the a-axis (see Fig. 2). In 1 
\ tr iilin 
\ P* —— ; 
ae ery | exper 
= 
\X 
° \ 
y*} A 
\ 
Ko > 
¥ sin ce oa ist tcinaia nit 
Le X " i a 
Fic. 2. Jet mean line geometry } 


The integrands of (8.1) and (8.2) may be written as (cos #) 4 and — (sin é 
respectively, where in the neighbourhood of P*, h = s = — R*(r*—7 


Thus the orders of infinity of the integrands of (8.1) and (8.2) are 


cos T* sin 7* 
sos cpenymenenm and paces . 
R*(7*—7) R* (7* — 7) 
and these functions are used to ‘subtract out’ the infinity of the integrands 
Hence, the value of the integrand of (8.1) at 7 7* is given by 
: X*—X cos T* the y 
lim > -\o se 7\> | Ds . lation 
roe | (X*#—X)?+(Y¥*—Y)?  R*(7*—7) oer 
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lues ¢ = y tT TE; 


| so that X = X*+€X_+4e2X,,+ O(e) 
| Y = Y*+c¥.+}4eY,,+O(e), 
s form | chere x. Reost, Y,= Rsinr, 
big X R_cost+Rsint,  Y,, = R,sin7-+Reosr, 
, the limit becomes equal to 
sinz* cos7dR 
5 WOl | 2R 2R? dr|\__. 


nteg? 8.1) may then be written in the form 





X*_X cose | 
. } \(x*#—XyP+() Y)2" R¥(r*—7)) R* 7 


in rly ~~?) may be written in the form 


Se i y*—) sin 7* lo < T, 


| (XPV RF) Re = 


), Numerical results 


In the rical example it was decided to take the jet near to the 
x edor order that the result obtained could be compared with 
| theoret results (4) already available at the time. 
as anes 
ee -—---- 
) Te et as 
; \ ration (To DDD, Cy 0-5 

55°5° and / 0-991 were therefore taken in our calcu- 


ond to an early model which had then been tested at the 
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N.G.T.E.; these experiments have now been superseded by more accurate 
tests on a closely similar model with 7, = 58-1° (7). 

Theoretical values were obtained for the lift, pressure distribution, and 
‘pressure-force’ contribution to the moment corresponding to a jet coeffi. 
cient C; = 0-5. The initial curve, marked 0 in Fig. 3, for the limit line of 
the jet as given by the corresponding source type flow taken from reference 
(6) was found to be a poor approximation to the true mean line, and the 
iteration was repeated four times, marked 1, 2, 3, 4 in Fig. 3, before satis. 
factory agreement between successive lift coefficients was obtained. Th 
final theoretical value of the lift coefficient was 2-60 for 7, = 55-5° as 
compared with an experimental value (7) of 2-9 for 7, = 58-1°, and with 
the value 2-75 for 7, = 60° obtained by the theory of reference (4). The 
pressure distribution around a 12} per cent. thick aerofoil was calculated 
using the usual transformations from flat plate to circle to ellipse, and is 
shown in Fig. 4. The pressure-force contribution to the moment coefficient 
about the leading edge was calculated to be 1-42 and the centre of tota 
lift at the point « = 0-55. 

_ 

— 

8-0} 
7-0} 
6-0} 
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Fic. 4. Pressure distribution on a 12}°, thick ellipse 
(To 55-5°, Cy 0-5) 


10. Conclusions 

An iterative method of solution for the equations of the limit line of th 
jet and consequently for the pressure distribution, lift, moment coefficient 
and centre of total lift, is given when the jet is at a general position on tht 
lower surface of an aerofoil. The length of computation involved make: 
the method too cumbersome for general use unless the initial curve chose! 


to represent the limit line of the jet is very much nearer the true positid! 
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1an that taken in the example described in sections 6-7, namely, the mean 


ne of the corresponding source type flow. A possible alternative would 


beto obtain the limit lines of the jet experimentally for a range of 7, for one 


ommon value of C,. For a given 7, limit lines of the jet for successive C,’s 


could then be obtained quite easily and quickly by iteration from the 


experimentally obtained limit lines. 


l 
2. F 
3, G. I. Taytor, Offerts ad Riabouchinsky. Memoires sur la Mecanique des Fluides. 
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LAMINAR BOUNDARY LAYERS AT THE INTERFACE { 
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SUMMARY 
The approximate solution of Keulegan (1) for the steady flow of a stream of viscous 
incompressible fluid over another at rest is extended to the case where both fluids 
are moving co-current but at different velocities. This solution utilizes a sext 
polynomial for the velocity distribution in the boundary layers. The solutions deper 


only on the ratio U,/U, of the velocities of the two streams and on the product of tl 





corresponding viscosity and density ratios. Numerical results are given for sever 
values of p12 p2/(p4, p,) at one value of U,/U,. Lock (2) has published an exact solutior 
with a numerical result for py po/ (M4 py) 1, U,/U, 0-501 and the sextice polynom 


solution is evaluated for this case also. Comparison indicates that in general the sext 


polynomial is more accurate than the quartic polynomial but that the advantage is 
not great. 


1. Introduction 


DISTILLATION, gas absorption, and liquid-liquid extraction, three important 
chemical engineering operations, are all instances of mass transfer from one | 
fluid phase to another, both phases being in motion. Little is known of the | 
fundamental hydrodynamics of fluid-fluid interaction even where stability 
of the interface is assumed. Keulegan (1), utilizing a sextic polynomial for 
the velocity distribution in the laminar boundary layers, gave a solution 


applicable only in the case where the lower stream was at rest. Late! 








Lock (2) obtained an exact solution applicable to the two cases—lowet 
stream at rest, lower stream in motion. Lock also gave the momentum 
integral solutions for various assumed forms of the velocity distribution 


including a quartic polynomial form which is the same as that used by 


Keulegan for his first approximation. 

These various solutions to the laminar boundary layers give rise t 
corresponding solutions for mass transfer between the two fluid phases 
Solutions to the mass transfer case have been obtained and will be pub- 
lished later. 

Unfortunately such solutions can be only of limited value since one of 
the prime factors in fluid-fluid interaction and in the behaviour of fluids 
with a free surface is the question of the stability of the interface or free 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 3 (1957)] 
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arface. Various aspects of this question have been attacked by Lessen (3) 


nd Yih (4), whose results confirm the expectation that instability sets in 


t rather low values of Reynolds numbers. 


2, The boundary layer equations 


(Consider two incompressible fluids flowing parallel and co-current; let 


the upper stream be of density p,, viscosity u,, and have velocity U,, and 


the lower stream have density py, viscosity py, and velocity U,. 


The origin is taken to be the point at which the two fluids first come in 


ntact. the axis of x is horizontal and in the direction of motion and the 


xis of y vertically upwards (see Fig. 1). The corresponding components 


hy \ 


| | 
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velocity are denoted by wand v. The boundary layer equations are then: 


: Ou cou C7u 
pper fluid u v v : 
as 
ca oy coy” 
° 
ou cu Cru 
wer fluid Vo , 
Ca cy ~ CY" 
here vy, and v, are the kinematic viscosities of the two fluids. 


The equation of continuity 1s 


0 
CXL cy 
+} ] ~ / 
iat there exists a stream Tunction w such that 
Cus Cus 


Introducing the dimensionless variables 





(1) 


bo 
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a solution is sought for in which: 


upper fluid yb = (v, U, x) f,(,), (6 
lower fluid yy = (v, U, x) fo(0). (7 
Then we have 
upper fluid u = U,f3(), (8 
lower fluid u = U,f3(ne), (9 
where accents denote derivatives: 
. , W/(Uyry\b¢ ps \ 
upper fluid oo N 2 tS 1(%1)—filai)}} (10 
<— : , 1/U,v2\3 i ’ ) 
lower fluid v= on {not o(n2)—fo(n2)}, (1 
4 7 Liat 
upper fluid i = l ( *) En) (12 
ey v4 
; eu + { U;\4 gn 
lower fluid — = U,|—)*f2(n-2) (13 
cy Vga 


Equations (1) and (2) then reduce to 
oh,» Mh, 
dy} ‘ ‘dn? 


9 Ufs Lf d*f 


and 2 f= 
dys * “dn 


= 0 (lower fluid). (15) 
2.1. Boundary conditions 
The boundary conditions (discussed in some detail by Lock (2)) which 
completely determine the solution of equations (14) and (15) are 


u>U, asyn,>, ie. fi(00) = 1, (16 
u—>U, as n> —0, i.e. f,(—oo) = U,/U, =); (17 
since the motion is steady, we have 
f,(0) = f.(0) = 0. (18) | 
For continuity of velocity, 
£1) = (0), as 
and for continuity of tangential stress, we have 
prrifi(0) = prvi fs(0). = 


2.2. Exact solution 


It is necessary to know the asymptotic forms which the solutions of the 


2f"(n)+F(a)f"(n) = 0 


equation 





take wl 


where 4 


where / 


with ec 

Lock 
and jy 
at rest) 


3, Apy 





= 0 (upper fluid) (14) | 
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take when 7 is large. For the case A > 0, Lock (2) wrote 


7 lim (f—An) E, 


shere E is a constant, and for 7 large and negative obtained 





se - 3 15 
fw Ant E+——{1-— 4 i— " (21) 
: 1 k2e-2=" 
f’ ~A+2kA erfz—s- —. = err (22) 
f” ~ kre-*—4kPat +... (23) 


ere k is a constant and 


rith corresponding expansions for f(7) when 7 is large and positive. 
Lock then obtained by numerical integration a solution for A = 0-501 
yay} ald pug Py/ (My Px) |. He also obtained solutions for A = 0 (lower liquid 


trest) at various values of py py/(u, p,), namely 5-965 x 104, 100, 10, 1. 





3, Approximate solution by the momentum equation 
i4)| Ifthe equations of motion (1) and (2) are integrated with respect to y 
etween the limits 0 and 6,, 0 and —6, respectively, where 5, and 6, are 
the boundary layer thicknesses on either side of the interface, then the 





1d ; ; 
mentum equations are obtained: 
01 
T) ou ¢ : on 
_ | upper fluid vs| -— | (U,—u)u dy, (25) 
) whicl Py OY) os64 Ox J 
0 
02 
16 : 2 2 To ou c U7. 94 
wer fluid vol (U,—u)u dy, (26) 
= Po “\ OY) y+ Ox . 
(ii 0 
vhere 7, is the skin friction. In the upper fluid we write 
Q u , 
] = $,(7}), 
1 
19 rere UB Y 0}. 
n equation (25) becomes 
1 
1 F dé < 
Z L d, 1) on 1 | (4, p*) dy}. 
l 1 dx J 
0 
1 
soit Pv. a ‘ 
ad j 9 » or 
ce 7 P,(V) OF | (4,—93) dnt. (27) 
1 R 


0 
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° UL 
For the lower fluid 7 = ho(7%), 
1 
where ns = —Yy/do, 
and equation (26) becomes 
| 
2v. XL ae a9 % 9 * 
TT $2(0) = 3 | (Ad.—¢3) dy. (28 
1 “ 
0 


3.1. Sextic and quartic polynomial velocity distributions 
The velocity distribution used by Keulegan (1) was a sextic polynomial 
which is again used here: 
; ~ _ l—dy * 
upper fluid ¢,(n*) = a+ - 9, (7%), (29 
a 


where, if the interfacial velocity is v9, ay) = Uy/U,, and also A, is a constant 


and 
* * ; *3 55 1 
91(ni) = (L+-my)ni—(1— Yl, + 10m, n+ (3-81, +9, ) nt 4 
+ (181, —27m,) nF? +9m, n**, (30 
where /, and m, are constants. 
Similarly for the lower fluid 
ioe A—a, 
$2(72) = A+ q: 
where A, is a constant and 
Jo(n3) = (1-+-mg)n3—(1— 3. + 10mg) n+ (—B8lo+ Ymg)nz + 


+ (151, —27m,)n*°+9m, n2*, (32 
where /, and m, are constants. 


If the constants / and m are put equal to zero then (as will be seen later|| 


A, } — A, and equations (29) and (31) become 
Pi(NT) = A+ (1—a)(2yf—2nP* + nF), (33 
$2(72) dy+(A Ay)(23 — 23° ne): (34) 


3.2. Boundary conditions 
The velocity distributions must satisfy the boundary conditions 
¢,(1) = 1, $,(1) = 4)(1) = ¢7(1) = 9, 
d(1) = A, $3(1) = 43(1) = 45 (1) = 0, 
$,(0) (0) = dp, 


6;(0) 3(0) = 0, 


Hay 
4 4’ (0) 
) 


_ He gr 0). 
0; 


> 
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Also, for »* = nf = 0 these differential equations must be satisfied 


when differentiated with respect to y. 





' 
3.3. Application of boundary conditions 
Since ¢,(1) = 1, then from equations (29) and (30) 
l—a —_ 
28 l a, °(4+81,+43m,). (35) 
~ 
Hence, A, = (1+ }81,+m)). (36) 
Similarly, fo1 the lower stream, since $,(1) A, 
_ A—a _ 
\ t °(4+81,+4m,). (37) 
ee 
9 Hence A, s(1-+ 18], +My). (38) 
Then, since p, (1) 0 = 5 (1), 
nstal . - 
t follows that 6 321, 60m, 0 (39) 
and 6+32/,+60m, = 0. (40) 
Using the condition for the continuity of stress, viz. 
Ni a Lo yr 
(0) Be 4° (0), 
O71 02 ~ 
u,(l-+-em | —@) t.(1+m,)(A—a 
ve have mM o) _ _ Bel 2)(A~4o) (41) 
0; A, 6, Ay 
If the equation of motion (1) is differentiated with respect to y and the 
onditions at the interface 1 n> 0 inserted, then 
‘i .- =H | dd l—a X 
Ly f °(1+-m,){ —8,— Vy; °(— 6+ 321, —60m,) 
A, ; dx A, 
nl : 
which, on integration, yields 
] a 2y, x(1 a 
Ao "(1+m, _ : 0) ( 6+-32/,—60m,). 
: A, 63 A, U, 
Using equation (39) this equation simplifies to 
641, a,(1- m,)dz(U, 2v; 5. (42) 
Similarly, for the lower stream, we have 
641, ay(1+-m,)d3(U,/2v, 2). (43) 


3.4. Integration of the momentum equations 
The sextic polynomial velocity distributions are used to integrate the 


omentum equations (27) and (28). For the upper stream 


9 


2v, x(1—a, ay(1—ay) Bb, o,(' iF 


. (44) 
67 U, A, A, aoe | 
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where B,(i+m,) = 0-15+-0-26667/, +-0-214285m, (45) 
and 
C,(1+m,) = 0-029365-+-0-07621, + 0-07626m, + 
+ 0-0309/? + 0-076191, m,+-0-0401 mj. (46 
Similarly, for the lower stream, 
2v_ x(A—4) =~ A—A , 10 (*- =] 1“ 
63 U, A, ; Ay : : A, , | 
where B,(1+m,) = 0-15+-0-26667/, + 0-214285m, (48 
and 
C,(1+m,) = 0-029365+-0-0762/, + 0-07626m,+ 
+ 0-0309/3 + 0-07619/, m,+0-0401 m5. (49) 


3.5. Method of solution 


For convenience we now introduce dimensionless boundary layer thick- 


nesses €, and €, such that 


») 2 } 
, [2,2 
1 
& 5,( q | 
1 
» . 4 
» [2ve2 
2 3, : 
UY 
l—a 
and put a, — 
A,é; 
ss A—ay 
r A, &, 
These definitions of a, and a, may be resta‘»d as 
ay+a,A,é, = 1, (50 
Ay +a, Ay& = A. (51 
Equations (41), (44), (47) then become, respectively, 
vos eo 
Ly Po\4 m 
a, (H2e2)*( ay aa (dla 
My pi) \l+m, 
2 242.1 (1 93 ¢3 
@; By ay ay ft +O, a; &, (448 
a} = Byaya3+Cy a8. (47 


Equations (42) and (43), together with (36) and (38), become respectively 


4a,(1—a,)? 


641, 5 [( . m,)(1+1381,-+-m,) 2), (42a 
ay 
4 —(§,)* , 
641, - apa [(1+m,)(1+1$1,+-m,)-*]. 43a 
as a8 2 
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(49 These equations, together with (39) and (40), constitute nine equations 
to determine Ap, 41, Hes S45 Sa: bs, my, l., and Mo. 

To obtain the first approximate solution the constants 1,, m,, l,, and mg, 
re put equal to zero. This gives the solution for the quartic polynomial 


velocity distribution—equations (33) and (34). The values of do, 4, A, E4, 


nd €, so obtained are then used to calculate values of 1,, m,, /,, and m, to 
‘7 eused in the second approximation. In evaluating /, and/, from equations 


) 


12a) and (43a) squares and cubes of the /’s and m’s can be neglected. 


3.6. Boundary layer properties 
The properties of interest are the interface velocity uy which is given 
lirectly by the solution of the equations (this quantity is highly important 
19 for the mass transfer solutions since it determines the time of contact), and 
he skin friction coefficient ¢ 


If the tangential stress at the interface is 7,, then 


ou a,(l+m,) of U, 2\-3 
Ty 1,\ -) 5 Mp, UI 7 > 


4, Results 
Calculations have been carried out for the cases 
A 0-501 and = py po/ (py py) a 
so for A = 10 and py ps/(, py) = 0-01, 0-10, 03162, 1-0, 3-162, 10, 100. 
In liquid-liquid interaction the values of py ps/(u, p,) of interest lie in the 
region about 1.) The results of these calculations are presented in Tables 
ind 2 
The first calculations made for this system are those of Keulegan (1) who 
resented the results of calculations for quartic and sextic polynomial 
city distributions for values of A = 0 and py po/(p, p,) 0-01, 0-10, 
03162. 1-0. 3-162. 10-0. 100. 
Lessen (3) appears to have obtained an exact solution for 
\ v M2 Po) (My Pi) l, 


» Lock (2) has caleulated exact solutions and approximate solutions 


juartic poly nomial or sine velocity distributions) for 
134 () Ls Po (fy P}) 5-965 « 104. 100. 10. 1. 


) ) | 
Po 1 Py : 
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TABLE 1 


Table of constants 1,, 1, and m,, m, 






























HeP2 

A | MaP1 I, My lg Mo 
= = | -_ = 

I'o O°1124 O°0401 0*1405 "0249 

oro! 1534 "O152 070568 0°0697 

or10 01661 O14 0°0794 0°0576 

Io°o 0° 3162 1707 0°0059 o*OgI2 O°0512 
I I°0o 1741 "0072 O'IOI6 | o70455 
I 37162 | 0°1763 0060 | o*10906 00416 
I Io°o 1777 0°0052 Ov1153 0255 
I T00°0 o*179¢ 0°00 O'r2I1 0°0354 









TABLE 2 


Boundary layer constants for the sextic polynomial solution 









HePe2 























A bP} 1g 1 as é, &5 
O'5o1 I'O 0*70600 Or1750 O1731 | 3°1005 #5555 
TO° Oo! 2° 299 1° 3865 14°634 | 2°2901 1°21 
: “¢ 2-889 -7903 | 12°573 | 18742 | 1°1348 
I 3162 | 4°9809 5°3197 10-812 | 1'69I0 | I0g05 
IC 1‘0O 6°1699 8°305 S-6411 1°5424 1°o481 
10°C 2162 72932 10°894 6°3541 1°4334 I°O123 
Ic 10° 8-217 12188 $°3150 1*3593 9955 
I I C 9° 3342 IO"1560 1°6672 1°2523 )501 

TABLE 3 
Interface velocities (u/U;) 
MeP2 Ilecurate Quartte Sextt 
A Py olution polynomial | polynomial 
5°965 104 O'O21I1 0°0220 
O 100 1727 01770 0°1739 
Q 10 3425 0° 3450 3440 
i I 5873 0*5907 g8ge 
o'sol I 7657 07668 7660 








TABLE 4 


l'angential stress coefficients (c,,) 






















I 





YI2190 


HePe Accurate Quartic | Sextice 
A Hyp solution | polynomial | polynomial 
° 5965 x 104 0°3318 | 0*3309 
fe 100 0°3183 0° 3186 0°3144 
fe) IO 02815 0°2774 2772 
° I o*1995 | 
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LAYERS 


In Tables 3 and 4 the various methods of solution are compared as far as 
existing calculations allow. It will be seen that the interface velocity given 
by the sextic polynomial is in every case better than that given by the 
juartic polynomial but this is not always the case for the skin friction 


oefficient c_.. 
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THE FLOW OF FLUID ALONG CYLINDERS 
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SUMMARY 

The boundary layer equations for uniform flow parallel to the generators of any 
cylinder without corners are put into the form of a series of linear third-order differen. 
tial equations. The first three of these are the same as those obtained by Seban and 
Bond (1) for a circular cylinder and solved by Kelly (2). The rest have additional 
terms depending on the radius of curvature of the cylinder and its derivatives. Th 
problem is also attacked by a Pohlhausen method as far as four terms of the series. For 
large distances from the front, Rayleigh’s method, as given by Hasimoto (3), gives 
the first two terms of an asymptotic expansion for the drag. Explicit calculations ar 
made of the drag of an elliptic cylinder of eccentricity v3. There is evidence that 


the drag is everywhere less than that of a circular cylinder of the same perimeter 


1. Introduction 
WE consider the flow of an incompressible fluid past a cylinder immersed 
in it, the fluid moving with velocity U parallel to the generators at infinity. 
Seban and Bond (1) considered the case of a circular cylinder and showed 
that, provided that the radius of the cylinder is large compared with the 
width of the boundary layer, the flow is the same as that for a flat plate; 
this is usually called Blasius flow. [See, for instance, Goldstein (4).] Since 
the boundary layer thickens as we go downstream, Blasius flow applies 
only for a short length of cylinder. By means of a suitable series, Seban 
and Bond extended the distance of applicability along the circular cylinder 
a similar attack on the problem has been made by Sowerby and Cooke (5). 
Glauert and Lighthill (6) obtained a Pohlhausen solution and also extended 
the distance still further downstream. They showed that the Pohlhausen 
solution is reasonably accurate. 

This paper treats the problem ot a general cylinder which has no sudden 
changes of curvature along its cross-section boundary. The problem of the 


boundary layer along a wedge remains unsolved. The recent attempt bj 


Loitsianskii and Bolshakov (7) only applies to the flow along the inside oi 


a wedge; the problem appears to be one of remarkable difficulty. 

We give the solution in the form of a series, the first three terms of whicl 
are the same as for a circular cylinder of radius equal to the local radius 
of curvature. It was, of course, to be expected that the first term would 
be the same (corresponding to the Blasius flow), and perhaps also the 
second. The succeeding terms consist of those for a circular cylinder plus 
certain extra terms which depend on the radius of curvature of the cylinder 


[Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 3 (1957)} 
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nd its derivatives. By the numerical solution of a series of third-order 
ferential equations general results in terms of certain parameters could 
egiven. Though the equations are linear, apart from the first one, the 
re olution is laborious and is not attempted here. 
However, a Pohlhausen method similar to the one employed by Glauert 
| Lighthill can be applied to the problem to give results in the form of 
series. It is not easy to gauge the accuracy of the results by this method, 
t the error may be expected to be much the same as that suggested by 
‘lauert and Lighthill in their solution of the problem of the flow along a 
Jar cylinder. Results are applied to an elliptic cylinder, and there 
ppears to be some evidence that the drag may be expected to be every- 
where less than that of a circular cylinder of the same perimeter. 
os The important results are that the first three terms of the Seban and Bond 
: series, and of the Glauert and Lighthill series, for the local skin friction 
neter fa circular cylinder of radius a apply to a general cylinder if the radius a 
sreplaced by the local radius of curvature p. Indeed, if the total drag per 
nit length of a closed general cylinder is required, it is possible to make 





rersec’ | use of four terms of the series for a circular cylinder. It must be emphasized, 
inity. \ however, that the general cylinder must have no sharp edges or sudden 
howed | changes in curvature. 
th t 
plate; | 2. The equations of motion 

Since The equation for the steady flow of a fluid is, in the usual notation, 
pplies q Xw = grad{(p/p,)+4q*}+vcurlw, 
Sebar 

dor Where w = curl q, q is the velocity and p, the density. 
a " We take an orthogonal coordinate system 2,y,0, where x is distance 
. ips | measured along the axis of the cylinder, the planes 2 = constant being 
oe | right cross-sections; @ is the angle which a plane through a generator 
ae | normal to the cylinder makes with some fixed plane belonging to the 
“e | family, and y is distance along a normal to the cylinder. The cylinder thus 
of | ngs to a family of parallel surfaces, y = constant and we may take 

om, 0 to be the given cylinder. If in the section by the plane = constant 
A the radius of curvature of the curve y = 0 is p, then in general we can take 

| the square of the line element to be 

whi | ds° dx*+-dy*+-(p+y)? dé. 
radius | We shall denote the three components of velocity corresponding to the 
would | three coordinates by u, v, w, and will assume that w = 0, so that the fluid 
30 the | tows along the cylinder without cross-flow. We shall also write r = p+-y; 
r plus | pis clearly a function of 6 only. 


linder | Wemay now derive the equations of motion, making the usual boundary 
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layer approximations, but assuming that r is the same order as 54, insteaj 
of order unity, as is usual. We shall also suppose that p’/p is of order unity 
dashes denoting differentiation with respect to 6. We then obtain 


ou, ou lop, vf{ée/ ou é [1 ou\) 
u—+tv—= —— —+-{_Ir—]-+ - —}}, 
Oxy py ex | rley\ ey) | e0\r e6)) 
0(8) + a 
Po CY 
l op 
O(o) = — —_ 
rpy OO 


Here 6 is the boundary layer thickness. 

The second equation shows in the usual way that the pressure may bx 
taken as constant across the boundary layer. We therefore take tl 
pressure in the boundary layer to be the same as in the main stream 
which in our case is constant. The third equation shows that in order for 
the flow to remain ‘straight’, i.e. w= 0, a cross pressure gradient i 
required. However, since the total change in pressure right round the 
cylinder is required to be of order 5°, our assumption w = 0 is justifie 
when the pressure is in fact uniform. 

Thus the equations to be solved, with uniform pressure, are 


ou, eu v ¢ ou\ v d/l du 
u—-+v— = - —|r—)+-—|- —], 
Cx cy r cy\ cy r c@\r 06 


= (ru) i. (rv) = 0, ; 
ox oy 


where r = p+y and p is a known function of @ only. 

For flow past a circular cylinder the equations are of the same fort 
except for the last term in equation (1). This term thus represents thi 
effect of the departure of the cylinder from the circular form. 


3. The solution of the equations 


We can satisfy equation (2) by writing 
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fre We write ‘ 4 (=) ; : l ("-)' (r2— p2), 
. ' 4p \vx 

b = p(vxU)'{ fo(n, )+EA (0, 9) +E*fo(0, 8)+.--}- 
[he velocity in the mainstream is U. 

We shall adopt as being better known the substitutions of Seban and 
Bond (1) rather than those of Sowerby and Cooke (5). They differ only 
y numerical constants and certain signs. We shall take y and (U/vz)# 
to be positive always; p and € will then be positive or negative according 

is greater or less than |p|, that is, according as the surface of the 
vlinder is convex or concave as seen from the fluid. The substitutions 

d to the equation 

u = WU (fottht+efet), 


lashes denoting partial differentiation with respect to 7. The boundary 





nditions are wu = v = 0 at 7 = 0, andu = U, v = 0 at » = ©, for all € 
mar nd 0. These he come 
ake 1 f (0 f.(0) 0 f(a) 2. f..4(®) 0 
SUrealll, | for g os. 2... . 
der f : 2 ‘ ‘ : 
‘ We express wu and v in powers of €, substitute in equations (1) and (2) and 
aient " > ° . ° . 
~ “t equate coefficients of powers of , making use of the expansions of p/r and 
und tl : 
pin powers of €y from the relation 
justine 
p (1- En); 
/ 
thisexpansion will be valid aslong as €7 is lessthan unity, that is, y << 0-414); 
| thus the boundary layer thickness should not exceed 0-414 times the radius 
| of curvature. With this proviso we obtain the following equations : 
lo JoJo Vv (3) 
| I Jo I lo hy 2hod Jo nto > (4) 
ae 2 Jole—=JoJot+%7o I nd fy T fy -2fi fi» (5) 
Ie it 
} " , L / ”" or gt or fn 
nts the} 43 +/ol3—Sfofst4fofs nds —fot+ file —3feh— 


—2fifet+A(n?fot4fi), (6) 








where A I (f _ 2? : L 
16 p p* 


Here primes attached to p represent differentiations with respect to 6, 






whilst primes attached to any f represent partial differentiations with 
respect to 7 
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Equation (5) should contain terms involving partial derivatives of / 
with respect to 0. However, since equation (3) does not contain @ explicith 


and since it must satisfy the boundary conditions for all 0, f, will be inde 


pendent of @, and so its partial derivatives with respect to that variabj 
do not appear in equation (5). From equations (4) and (5) the function; 
f, and f, are in the same way successively seen to be dependent only on 
Equation (6) will contain partial derivatives with respect to @ of earlie 
functions which also vanish. Higher order equations will be considerabjy 
more complicated, as they will involve non-vanishing partial derivative 
with respect to @ of some earlier functions. 

It will be seen that the first three of these equations are the same « 
those found by Seban and Bond for a circular cylinder. The fourth an 
later equations are the circular cylinder equations with extra terms ¢ 
pending on p and its derivatives on the right-hand side. A univers 
solution of ¢quation (6) is g3(m)+ Ahs(7) where 


M(D)gs = —nfe —f2+ 3h fa— fof —hfe 
A D)hs = Ai t+4fi 
MD) = PLP D+if, DP = o/0y. 


Each solution is to satisfy the same boundary conditions as fs. Similar 
methods may be used for the higher equations. 

Because of the difficulty of solving the equations of Seban and Bond }y 
desk machines, no attempt will be made here to solve them. Sower)) 
and Cooke found some instability in the solution of equation (4) by step: 
by-step methods. They obtained 0-692 for the value of f{(0); Seban an 
Bond gave 0-704, and Kelly (2), using an analogue type of computer, gav' 
0-696, which we shall adopt. We see that there is some doubt about t! 
correctness of even the second significant figure. Greater variations occ! 
in the case of /5. 

From the results we may find the skin friction 7 given by 

T pu(ou ey), 03 


we obtain 


T l ow i vr iad vies 
a = fo (9) +f (9) + €fo(0)+ €*193(0)+ Ah;(0)} 1 
“as S 


1328-1 + 0-696 — 0-199 + €{95(0)+ Ah3(0)}+-... 
using Kelly’s results (2), as far as they were calculated. 


Because of the difficulties mentioned relating to the solution, we shi 
attack the problem by Pohlhausen’s method. 
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.[4, A Pohlhausen solution 


| Following the method of Glauert and Lighthill (6) we take 


U \ 
u log(1 m 5 
Xx \ P) 


sour distribution, where « is now a function of x and 6. 


(3) 


This form satisfies equation (1) at the wall of the cylinder, where y = 0, 
but not the equation obtained from equation (1) by differentiating it with 


Irespect to y and putting y = 0. It can be shown, however, that the un- 


satisfied terms will only be significant if the series which we shall obtain 
iter for € (see p. 319) is carried to higher powers than is done in this paper. 
The edge of the boundary layer is given by y = 6, u = U, where 34 is 
Hence we have 


x), 


the boundary layer thickness. 


It will be noticed that this velocity distribution does not satisfy the 


ndition éu/éy = 0 when y = 6. Glauert and Lighthill consider it impor- 





tant to fulfil as many boundary conditions near the wall as possible, and 
Itheir results justify this conjecture. 


We next obtain a momentum integral equation. We do this by multi- 

} ; : 
plying equation (1) by r and integrating with respect to y between the 
) 


inits 0 and 6, making use of equation (2). The process is well known and 


‘given in Goldstein’s book (4). 


The result is 


0 


: eu) U. 
: u(U uy)? dy oF, me. E (9) 
CL. eyo x 

; ts) 
F - C 
where 1 : ( aa dy. 
¥ J fale] \ 7 fale] ; 
0 


| The first term in equation (9) has been shown by Glauert and Lighthill 


0 to be 
Dw2 2 9 2a 412) Z 
| prom ¥ mdb SS. (10) 


4a3 Ca 


the term —v{r éu/éy}2 is taken to vanish at the upper limit (as should 


de the case) though in fact it does not do so because, as has already been 


pomted out, the particular boundary condition @u/éy = 0 when y = 8 is 
not satisfied by our form (8). At the lower limit the term becomes vU/a. 
The evaluation of J is straightforward. We have 
C {1 eu\ 
rie aad 


“a a 7 


f- 
<— + (Bp?—B’) . log( 2. “| 4 “Flog (1 a's *). 
; p) p 
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where a = p’/p, 8 = «’/ax, dashes denoting derivatives with respect to 
and keeping x constant in the case of « and f. 

Integrating this expression between the limits 0 and 6, that is, betweq 
0 and p(e*—1), we have 


| YT dy = —1+a+e-, | od dy = 3(1—e-*)?, 
2 ;: 


[ “log{ 1 +. ¥) dy = }a?, | P log(1 + A] dy 1—e-*— ne-*, 
r p = a oa . p 


and so 


J = aB(—1+2a+e-%— ae“) + a?(1 —e-*)?+- 407(B?—B’)—a'(—1+at+e 


« 


From equations (9) and (10), on multiplying by 4a/p?U?, we have 


D9 y2—§ ah.  -— 19) ¢ . 
te ee. #4 r 


ss Cx p* 
where J is given by equation (11). 

The solution of equation (12) is in general complicated, but it can | 
obtained as a series.t We note from equation (11) that J is of order « 
Hence we first put J = 0 in equation (12), obtaining as the first approx 
mation the result of Glauert and Lighthill, viz. 

Dur! n2y7 2 4.3 
€ (= 12va/p?U) 05 —f-$00* + 04, (1 
the series being correct as far as written if J is of order «?. 

On differentiating the series with respect to @ and rearranging, ¥ 


obtain , In 2 
Cu aaa“ 


oo 3 
and since B = a’/a, 


, , , 9 
B a+-—(a—a*)+ 


Substituting these values in the expression (11) for J we find that t! 
coefficient of «? in J vanishes and that the coefficient of a* is 4a°—; 
Consequently, as far as the term in a’, the value of J is —8Aa°/3 in‘ 
previous notation. Hence J is of order a® and so we may write anotli 
term of the Glauert and Lighthill solution into the series for «. Henet 


2.14.3.) 9.4.1 51 14 
€ x Bie 10% As a eeey 


+ I am indebted to a referee for simplifying the original procedure adopted in solving! 
equation. 
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ct to (} where a, is a function of @ to be determined from equation (12). Differen- 
tiating the series (14) with respect to x, using equation (13) and substituting 
et weeny in equation (12), we obtain 
Mw, 3, 8A > ag ae P 
4o°+-— a —-a*+... 14 x? +-...]{ 2a+ 4a°+ — a? + 5a, a*+...]. 
» 15 3 v 


The coefficients of «, «?, and a® agree, as was to be expected, and we 
1y determine a, by equating coefficients of a4. We obtain 


) 


| o - ) : 
€ x" + . ?. x + — x” “T" ese 


715 15 
This equation can be inverted to give 


l 2 13 fe is) 
j by . e- 


, 3. «660 675 15 


using our previous notation, 


1-188€-1+- 0-667 — 0-188 + (0-149—0-400A )é?+... 





rder } which should be compared with equation (7). Glauert and Lighthill gave 
pp! series similar to (15) as far as the term in e!. Their next term would have 
been 134€/675, that is, it would not have included the part —8Ae«/15. Once 
more we note that the first three terms are the same as for a circular cylinder, 


ut that there are differences in the fourth and succeeding terms. 
5. Drag for small a 
Since the drag per unit area is r, the drag D per unit length is given by 
D | 7p dé, 


the limits of integration being 0 and 27 for the whole of the outside of a 


osed cylinder. In general the drag is thus given by 





; D = pl § FZ, (16) 
n r=0 
where i? 6va/Ua?, a is a representative length, and 
*fq\r 
J | f” (0) dé. 
J \p) 
| We may note that since fj = g3+ Ahi, and 


9 


A 1 6p? “a(-) 
ds? p 
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where ds (= pd@) is the element of arc of the section, the value of J; for 
a closed cylinder reduces to 


| (2) 95(0) dé, 
lng 


since the integral of the term involving AA3(0) vanishes. 


6. Drag for large x 

It is not so easy to find a satisfactory method of dealing with this region 
but we may note that Glauert and Lighthill (6) and Stewartson (8) have 
obtained (consistent) asymptotic solutions of the problem of the circular 
cylinder, and that these agree in their first and second terms with the 
drag given by the well-known Rayleigh method in which an infinite cylinder 
is impulsively started from rest with a velocity U parallel to its length, 
Batchelor (9) made the first attempt to solve the Rayleigh problem for the 
general cylinder and Hasimoto (3, 10) considerably improved on this. W: 
therefore take the first two terms of Hasimoto’s solution in the expectation 
that they may be applied to give a correct result for the drag in our boundary 
layer problem as far as they go. 

Hasimoto effectively transforms conformally the section of the general 





cylinder into a circle whose radius c is chosen in such a way that the tw 














regions at infinity are unchanged. Batchelor calls the cylinder having this | 
section the ‘equivalent circular cylinder’ and he gives the first term for the 
drag. In fact, from Hasimoto’s work we now know that he could have 
taken the first two terms of the well-known formula for the drag of a| 
circular cylinder (Jaeger (11)) and written down the result for the drag D per 


unit length in the form D b + 


4mns:U BPP 
where B = log(4vt/c?)—2y, y is Euler’s constant, and c is the radius of the 
equivalent circular cylinder. 
Alternatively, keeping two terms in 5, we may write 


_D = a ae (17 
dal 8° 6 
where 5 = log(4vt/c?). 


Finally, then, on the assumption that Rayleigh’s problem may pass ove! 
into the boundary layer problem by writing ¢ = 2/U, equation (17) holds 
in the latter problem when we write 


8 = log(4va/Uc?). (18 
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7, Drag of an elliptic cylinder 
We will illustrate these results by finding expressions for the drag an of 
lliptic cylinder of perimeter 27a and eccentricity v3/2, that is, having a 

thickness ratio of 1/2. 

The numerical work involves the complete elliptic integrals EZ and K. 
The result is that 

D = 2npU{k- fo(0)+f7(0)+ 1-629kf3(0)+-3-576k293(0)+-...}, 
where k2 = 16va/Ua?, as compared with the drag of a circular cylinder of 
the same perimeter which is 
D = 2npU{k-* fo(0)+f7(0)+kfe(0) + k*g3(0)+-...}. 

It is to be understood, of course, that these formulae apply for small 
values of x. 

The first two terms are the same in the two expansions, but in view of 
the fact that f3(0) is negative the drag for small values of k is less for the 
elliptic cylinder than the circular cylinder. The matter is not so definite 
it larger values of k since g3(0) appears to be positive from the Pohlhausen 
solution given in section 4. 

For large values of x, if the cylinder has major and minor semi-axes a’ 
nd b’, the equivalent circular cylinder has a radius }(a’+-b’), whilst the 
perimeter of the ellipse is 4Ha’, which must be put equal to 27a. 

Hence the radius c of the equivalent circular cylinder is given by 
am{1+-(1—e?)!} 


 4éB 
where ¢ is the eccentricity and EF has modulus e. 
When ¢ v3/2 this gives c = 0-803a and consequently from equations 
7) and (18) the drag of the elliptic cylinder is again less than that of the 
ircular cylinder with the same perimeter. This can be shown to be true 
ior all elliptic cylinders, and indeed it seems very likely that the result is 
true for all closed cylinders. 
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SUMMARY 

The theory of simple harmonic electromagnetic waves in a nearly periodic structu 
of the type used in linear accelerators is discussed by an expansion in a series of 
appropriately defined transmission modes of the field in a unit cell of the structur 
The possibility of the expansion being assumed, it is shown that the coefficients car 
be expressed as integrals of the field over the input or output apertures of the unit 
cell. The modification to the corresponding solution in the strictly periodic structur 
can be calculated by perturbation theory. In the first approximation the structur 
can be represented by a series of four terminal networks, and the voltages and cur 
rents in these can be defined in such a way that the relevant parameters var 
smoothly in the neighbourhood of resonance. 
1. Introduction 
THE problem discussed in this paper is that of the calculation of simp 
harmonic electromagnetic fields in structures which are nearly periodic i 


the following sense. The type of periodic structure concerned consists of a1 


infinite series of identical cavities, each with perfectly conducting walk 
and an input and an output aperture, connected together so that the output 


aperture of each is the input aperture of the next. The nearly periodi 
structure is an approximation to a finite length of such a structure in whic 
adjacent cavities are nearly, but not quite, identical; the difference betwee! 


cavities separated by a large number of others may be large. The apertures 


of all the cavities, however, must be exactly congruent. 


Consideration of this problem is suggested by the use of such structure 


in linear accelerators for charged particles, which has been extensivel 
developed in recent years (see, for example, Slater (1), Alvarez et al. (2 
The present investigation was inspired, in particular, by the accelerata 
developed by Alvarez et al. (2), but is not confined to this case. The results 
are also applicable to an electrical transmission system consisting of a serie: 
of 2n-terminal networks, half the terminals of which can be considered & 
input terminals and half as output terminals. Such a network is, in fact 
with suitable geometrical design and disposition of the circuit elements, 
particular case of a structure of the type considered, the operating frequen¢ 
being low enough for the equations of alternating current theory to provid 
an adequate description of the electromagnetic behaviour. 

The characteristic property of a periodic structure of the type considered 


Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 3 (1957)] 
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sof all periodic oscillating structures (Brillouin (3)), is their property of 
sustaining oscillations in which the oscillatory field has the same form, but 
possibly different amplitude and phase, in all the cavities, the amplitude 
tio and phase difference between all pairs of consecutive cavities being 
Each 


transmission mode, considered from the point of view of a single cavity of 


the same. Such oscillations will be called transmission modes. 
the structure, is an electromagnetic field which has the same boundary 
values for its tangential components at the output aperture as at the input 
iperture, except for a difference of amplitude and phase. Thus one can 
speak of a transmission mode of a cavity as well as of a periodic structure 
provided that input and output apertures are congruent. (It is this which 
ecessitates the condition of exact congruence of all the apertures in the 
nearly periodic structures under discussion.) 

In the designs which have been used in practice the operating conditions 
we such that only one mode, or more precisely one pair of reciprocal modes 
nthe sense defined in sections 2 and 3, is propagated without attenuation. 
In a strictly periodic structure each of these two modes would have the 
same amplitude in all cavities and the same phase difference between any 
pair of successive cavities. In a nearly periodic structure in favourable 
circumstances only such dominant modes will be effectively excited, but it 
will still be necessary to calculate their amplitude and phase in the different 
cavities as in a structure with many cavities these may differ significantly 
from the unperturbed values. 

This calculation may be considered as part of a first-order perturbation 
theory approximation to an exact calculation of the field in the structure. 
Estimates of this approximate solution for the Alvarez accelerator have 
been made by a precisely formulated analogy with the effect of a small 

riation of diameter of an empty cylindrical waveguide (2). The problem 
has also been discussed by an analogy with one-dimensional wave propaga- 
tion in a continuous medium with slowly varying propagation constant 
Wilkins (4 

he method adopted in the present paper is an expansion of the field in a 
ity, or rather of the tangential components of this field in the input and 
utput apertures, in a series of transmission mode fields. The possibility of 


such an expansion being assumed, it is shown that the coefficients of the 


expansion at the output aperture for a given field there can be expressed as 
ntegrals over the output aperture, and that the coefficients at the input 
perture can Hence the field at the 


Nput aperture can be determined; and since this input aperture is the out- 


he easily determined from them. 


wt aperture of the preceding cavity, by repetition the fields at the first 


nput 


iperture can be determined from those at the last output aperture. 
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The expansion in modes is valid only if the aperture components 
of all the modes together form a closed set for electromagnetic field 
in the aperture. This is known to be so in certain cases, and it js 
assumed true for the cavities considered. The proof of the existence jp 
general of the reciprocal modes which are used in determining the ¢o. 
efficients in the expansion also depends on this assumption of closure, 
But in the case of unattenuated modes in all cavities, and all modes ip 
symmetrical cavities (except for certain critical cases), the existence of the 
reciprocal mode can be shown directly. When the reciprocal mode does exist 
the first-order perturbation theory approximation to the solution mentioned 
above is valid independently of the assumption of closure and of the 
validity of the expansion. 

By an appropriate transformation the relations between input and 
output fields for a single cavity can be expressed in terms of conventionally, 
but precisely, defined mode voltages and currents, so that each cavity is 
equivalent to an infinite set of independent four-terminal networks. In 
this notation, in the first approximation solution, the structure is repre- 
sented by a single series of four terminal networks, one for each cavity and 
one for each junction between cavities. For asymmetrical cavities the 
voltage-current representation is rather complicated, and the mode 
amplitude representation is better. For symmetrical cavities the voltage- 
current relations become simpler, and the junction networks reduce to pure 
transformers. 

The voltage-current representation is most valuable for the critical case 


in a symmetrical cavity, in which the dominant mode in the cavity is a | 


resonant mode with zero tangential electric, or magnetic, field in the 
apertures. The mode of operation of the Alvarez accelerator is of this 
type. In this case the dominant pair of reciprocal modes coalesces in the 
single resonant mode, and, since the tangential electric or magnetic field 
in the aperture is zero, the aperture fields of the modes no longer form a 
closed set. This necessitates a modification in the mode expansion to replace 
the missing term, which results in a discontinuity in the mode amplitude 
representation at resonance. The voltage-current representation, however, 
can be defined in such a way that all the relevant parameters vary smoothly 
in the neighbourhood of resonance, so that calculations are made mucheasiet. 

The present paper is confined to investigation and development of the 
general method and is of quite general application. It seems, however, t0 
provide a practicable method for computation, particularly of resonant 
structures such as the Alvarez resonator. An example of application to 4 
practical structure is given in section 7 to illustrate this, and further 
investigations are proceeding. 
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2. Transmission modes in a cavity 

Consider a cavity which consists of a finite empty region D, bounded 
partly by a surface S’ occupied by perfectly conducting material and partly 
by two apertures S, and S,, the input and output apertures. The positive 
normal to the bounding surface is taken directed outwards from D on S’ 
ind S,, and inwards to D on S,, and the apertures are to be congruent in the 
sense that there is a rigid displacement P — P’ which carries S, into S, and 
preserves the sense of the positive normal. In virtue of this congruence the 
cavity can form one of an infinite series of congruent cavities, coupled 
together so that the output aperture of each is the input aperture of the 
next, the whole forming a periodic structure such that the displacement 
P-» P’ carries each cavity into the position occupied by the next. 

The components Ee'¢*’, He’*“ of a simple harmonic electromagnetic field 
f wave number k/2z7 in D satisfy Maxwell’s equations 


VAE ikZ,H, VAH ikY, E, (2.1) 
in D, and the boundary conditions 
nAE 0 on VS, (2.2) 


where n is the positive normal unit vector on the boundary of the cavity 
ind Z, | Y, is the characteristic impedance of free space. (Rationalized 
m.k.s units are used.) 

If (E, H) is an electromagnetic field which satisfies (2.1) in the interior 
of the extended structure formed by the repetition of the cavity, and also 
satisfies (2.2) on its boundaries, then the field (E’, H’) obtained by subject- 
ing (E, H) to the displacement P + P’ is another. The whole set of such 
fields for a fixed value of k constitutes a linear vector space, and the equation 

T(E, H) = (E’, H’), (2.3) 
where E’( P’) E(P), H’(P’) H(P), defines a linear transformation, T, 
of this vector space into itself. An eigenvector of this transformation with 
eigenvalue m, that is, a field (E, H) such that 

T(E, H) (mE, mH), (2.4) 


is a field which is the same in all the cavities of the extended structure 
except for multiplication by a factor 1/m in going from one cavity to the 
next. The corresponding field in the cavity D thus satisfies the further 


bounda Vy ondit ions 


E(P) mE(P’), H(P) = mH(P’), (: 


bo 
or 
— 


where P is any point on S, and P’ is the corresponding point on S,. 
A field which satisfies (2.1), (2.2), and (2.5) will be called a transmission 
mode for the cavity, and m will be called the transmission factor of the mode. 
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In a vector space of finite dimension number every linear operator has at 
least one eigenvector. A 2n-terminal network with n input and n output 
terminals is a system with 2(n— 1) degrees of freedom, since the independent 
variables can be taken as the n— 1 independent input and output currents, 
Such a network therefore always has at least one transmission mode. | 
the dimension number is infinite, as is generally the case in the present 
application, there may be no eigenvectors, and it is therefore not obvious 
that all cavities of the type under consideration have transmission modes 
However, such modes are known to exist in many cases of practicai impor. 
tance; an example additional to those of the linear accelerators mentioned 
in the introduction is that of a finite length of a cylindrical waveguide ter- 
minated by plane apertures perpendicular to the generators of the cylinder, 
The well known waveguide modes are transmission modes for such a cavity 

3esides the transmission modes there may, according to the general 
theory of linear operators, be associated with a given eigenvalue, m, of Ta 


set of fields (E,, H,), (E., H,),..., (E,,, H,,) such that 


T(E,, H,) (mE,, mH,) | 
~ (am 
T(E,., H,) (mE, } E, ‘5 mH, H., 1) for +r 7 ie: SOE n | 
Such a set of associated solutions, with m +1,” 2, is of importance in 


the case of a symmetrical cavity operated at a resonance frequency. This 
case is discussed in section 6. 

To every solution (E, H) of (2.1) there corresponds the conjugate solution 
(E*, —H*), (* denotes conjugate complex). If (E, H) is a transmission 
mode with transmission factor m, the conjugate solution is a transmission 
mode with factor m*. If |m 1 then m* 1/m, and the conjugate 
mode may be considered as a mode with factor m in the sense from S, to 
S,, a mode ‘reciprocal’ to the original one. The existence of a reciprocal 
mode corresponding to any given mode in a symmetrical cavity (except 
for modes with m +1 in which the two modes in opposite sense may 
coincide) is obvious from symmetry. There is no such direct proof of the 
existence of reciprocal modes in asymmetrical cavities. It will be shown 


in section 3 to follow from the assumption of closure. 


3. Orthogonality relations between transmission modes 
If (E,, H,) and (E,, H,) are any two electromagnetic fields in D it follows 
from (2.1) that V.(E, \H,—E, \ H,) 0. If the fields also satisfy the 


boundary condition (2.2) the use of the divergence theorem gives the equa: 


tion 


( (E,.H,—E,.A,) dS = [ (E,.A,—E,.fi,) ds (3.1 


So Sp 
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here H = Hn. The integrals depend only on the tangential components 
f the fields in the aperture. If (E,, H,) and (E,, H,) are the fields of 
transmission modes with factors m, and m, respectively, (3.1) reduces to 


(I—m,m,) | (E,.H,—E,.H,) dS = 0, (3.2) 


where S denotes either of the apertures S, or S,. The surface integral is 
therefore zero unless m, 1/m,. 

Now consider the Hilbert space whose elements are pairs (A, B) of tan- 
sential vector functions defined on an aperture S congruent to S,, with the 


w product definition 


((A,, B,), (As, B,)} = [ (A,. A+B, .B%) dS, (3.3) 


S 


which A,. A denotes the ordinary scalar product of cartesian vectors. 


3.2) shows that the pairs of fields which consist of the tangential com- 


ponents in the aperture of (E,, H,) and (A%, —E%) are mutually orthogonal 
elements of this vector space if m, # 1/m,. The result obviously holds for 
ny value of m, if E, E, and H, H,. 
We now assume that the aperture tangential components of the fields 
fall transmission modes of a given wave number in the cavity form a 
plete set of elements of the Hilbert space; that is, that any pair of vector 
fields which is orthogonal to all of them is zero. This assumption is equiva- 
nt to the assumption of closure in the mean square sense. The investiga- 


tion of the validity of this assumption for a given cavity would be difficult. 
It is known to be true for cylindrical waveguides (subject to restrictions 
on the smoothness of the boundary of the cross section) except at the critical 
frequencies for the modes, and it seems not unreasonable to assume that it 


Wi 


| hold for any cavity of reasonably smooth shape except in certain 
ritical conditions of the type discussed in section 6. 

\pplying the completeness assumption to the element (H*, —E*) where 

E and H are the aperture tangential components of a transmission mode 

with factor m, we see that there must be a mode with aperture components 

t orthogonal to this element, that is, a mode with factor 1/m; this is the 

reciprocal mode. The two modes must be distinct even if m +-1, so that 
l/m, since (H*, —E*) is orthogonal to (E, Al). 

The transmission modes therefore occur in general in groups of four, 

each one being associated with conjugate, reciprocal, and conjugate 

reciprocal modes with factors m*, 1/m, and 1/m* respectively. The four 


may reduce to two if m is real (m m*) or unimodular (m 1/m*) since 
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in those cases a mode may be identical with its conjugate or with its con- 
jugate reciprocal. It-has been stated (Slater (1)) that in the case of electro. 
magnetic waves in lossless cavities m is always real or unimodular, but an 
example given in section 5 shows that this is not true. 

The transmission modes can thus be grouped into pairs of reciprocal 
modes. The pairs will be numbered 1, 2,..., and the members of the rth 
pair labelled r+ and r— (the transmission factors being m, and 1/m, 
respectively). For each mode choose a definite amplitude as unit amplitude 
and denote the tangential electric and magnetic fields in the output aper- 
ture for this amplitude by & and # with the appropriate suffix. The unit 


amplitudes of conjugate modes are to be chosen so that &F, = &,, 
H*, = —H,, where r’+ and r’— are the labels of the modes conjugate 


to r+ and r— respectively. Then by (3.2) the two sets of elements 
7 ae RN (6,_, %,_),... and (#*%_, —&*_), (#*,, —E&*,),... in that order 
are a pair of complete bi-orthogonal systems with respect to the scalar 
product (3.3). We define 

[ (6,,.%_—#,,.6,) dS = N, 


7 





 (6,...%,—#,,.6,,) dS = 0 : (3.4 
§ 

[ (6,,.%_-#.,.6,)dS =0 for r#s 

5 


Since the systems are complete N, is not zero; the unit amplitudes can be 


chosen so as to make N, = 1, but this normalization is not always con- 
venient. 
In view of the orthogonality relations any pair (E, H) of tangential 


vector fields in S can be expanded formally in a series of the form 


E = 9 (a,,¢6,,+4,_6,_) \ 
. : ; (3.5) 
A = 5 (a,, %,.+4 #,.) | 
ae ~ 
where a,, = + y (E.4%,-—H.4,-) dS. (3.6) 


It has been assumed for the sake of simplicity that there is only one 
linearly independent transmission mode for each transmission factor ™, 
(and only one pair of reciprocal modes for m, = +1 if these values occur 
at all). If there is more than one it can be shown without difficulty that 
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,linear basis can be chosen for the modes associated with each m, in such 
, way that the orthogonality relations (3.4) still hold. 

We now consider the expansion in a series of the type (3.5) of the aperture 
components of an arbitrary electromagnetic field (E, H) in the cavity, and 
in particular the relations between the amplitude coefficients a,,, and 
i,,, for the input and output aperture fields. The series (3.5) may in general 
have only formal significance, but under sufficiently restrictive conditions 
n E and H it will be convergent and will determine also an expansion 
f the field everywhere in the cavity as the sum of a series of transmission 
modes. The coefficients in this series will be the same as the coefficients 
in the series (3.5) for the aperture field in the output aperture S,, and by 


the defining property of transmission modes it then follows that the 
amplitudes for the field at the input aperture will be given by 


Bosc M, Ap. bs a A, p/| My (3.7) 


Relations (3.7) can, however, be obtained directly for the amplitude 
oeflicients as defined by (3.6) without reference to the expansion, and even 
without the completeness assumption, provided that the reciprocal mode 
exists and the normalizing integral N, is not zero. In fact equation (3.1) 
with (E,, H,) as the given field (E, H) and (E,, H,) as the r+ mode of unit 
amplitude, with aperture components (6,,, %,,) on S, and (m,6,,,m,%,.) 


rt r+ 


nS, gives immediately 


m, | (E.#,,-H.6,,) d8 = { (E.%,-f—1.6,,) dS, 
‘i Sp 


r from (3.6) N,Q, q = Ap—ps 


which is the second of (3.7). The first of (3.7) follows similarly by taking 
E,, H,) as the r— mode of unit amplitude. 

The amplitudes a,, and a,_ of the aperture components of a pair of 
reciprocal modes can be expressed by means of a linear transformation in 
terms of conventionally defined voltages v, and currents 7,, and the input- 

itput relations then assume one of the standard forms for a four terminal 
network. The definitions are 


v Bn Cys +p Cyns ip = @,,h,,+2,_h,_, (3.8) 
Where ¢,,., ,. are arbitrary constants, which for later convenience are so 
hosen that « ca h* , which may be interpreted as the 


oltages and currents in the rth mode for unit amplitude of the r+- modes 
ithe aperture. (3.8) can be solved for the amplitudes to give 


(3.9) 
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where Z,, = e,,/h,, are the characteristic impedances of the r+ and r- 
modes. The transfer relations for voltages and currents can be obtained 


from (3.9), (3.8), and (3.7) in the form 


m, Z,,—Z,_/m, Z,,4,(m,—1/m,) . 
v = —_—_—— v.. 5 = -v.} 
r.a 7, ,, r.b z, 7, r.b 
Z,+—4%, Ly +—Z, | 
. : (3.10 
, m,—1/m, m, Z,.—Z,,/m, . 
i ey 
r,a r ,, rb ,, 7 r,b 
Z,.—4, Z,,—Z, 
The constants e,,, 4, may be so chosen that Z,., Z,. = Z,, andif 


this is done (3.10) simplify to (3.19) below. 

The power transmitted through the cavity by a given electromagneti 
field can also be expressed simply in terms of the amplitudes or the voltages 
and currents. The mean rate of energy transfer P across the surface § i 
the positive sense by an electromagnetic field (E, H) is given by 

P=—1/| (E.A*+E*.A) dS. (3.11 
S 

Since the field conjugate to a transmission mode is also a transmissior 
mode the aperture components of the field (E*, H*) can be expresse 
explicitly as a series of transmission modes in the form 


YY sk . * 4 | * C 
E > (a? é,,1+-af_é,_), 


r 





A* > (at, #+0%_#_). 


Tr 


Thus, by the orthogonality relations (3.4), (3.11) gives 


L dw rp 
Po 


P a. N(a* a ax @,..,)- (3.12 


The various arbitrary constants have been chosen so that 


6. é* P A H* * ey e* , he, —h* 
from which follow (3,13 
Z, —Z* , N, = —N* 


[t can be verified that these conditions are compatible for modes which ar 
self-conjugate or conjugate to their reciprocals. With them and (3.! 
(3.12) ean be transformed to 
r L > {N(0, iF + v8 t,)/(e,.h,-—e,-h,)}, 
m 
so that if the constants are chosen so that 


N. e.,h,_—e,_h,.., (3.14 


r r+ r r r 


a choice which is compatible with (3.13), the expression for the pow! 
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nsfer reduces to P 


LY (v,8+043,), 


For modes which are self-conjugate or conjugate to their reciprocals 
,¢ i and the corresponding term is real and occurs only once 


the summation. Every other term appears together with a comple- 


entary term v, 7* +-v*i, which is conjugate to it. Thus the formula for 


ver transfer can finally be written as 


> 1 "a me oo hee 21K 
i i> v4, Ur ty) sre; > vir}, (3.15) 
r 
generalization of the formula P = }refvi*} for power input at a pair of 


rminals. 
[he introduction of the voltage-current notation does not simplify the 
wnsion formulae in the general case. For a symmetrical cavity, 
wever, i.e. a cavity which has a plane of symmetry such that S, and 
we images of each other in it in such a way that any point P’ of S, 
the image of the corresponding point P of S,, such a simplification is 
xsible. The simplification arises from the fact that to any solution 
E. B., H,, H,, H, of Maxwell's equations there corresponds a ‘reflected 
ition’, E’., E,, —E., —H H’,, H., where E’(x,y,z) = E,(x,y, —z) 
Ina symmetrical cavity in which the plane of symmetry is taken as the 
une : 0, the reflected solution satisfies the boundary conditions (2.2) 
he original solution does, and also satisfies (2.5) with constant 1/m if the 


ial solut 


tion does so with constant m. Hence the reflected solution 


transmission mode c: ve taken as the reciprocal mode. (This is 


r) 


l 


ious if there is only one mode for a given value of m. It can be arranged 


suitable choice of the linear basis even if there is more than one.) 
The fields and arbitrary constants of the reciprocal modes can therefore 
hosen so that 
( é Hf H i e.. é,, a 
} Rh. Z, —Z, Z, 
orthogon lity relations (3.4) reduce to 
— o£ ra, ea 
6..%. 48 ‘ a (3.16) 
, | LN eh. if s - 
\ ~ / 7 7 
here e. and h. have been chosen so that N, has the value (3.14). 
The expansion relations (3.5) and (3.6) reduce to 
; a ~ a — 
g- 574, 8a 5s. (3.17) 
iat ah. 


r 
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and the transfer relations (3.10) to 


_— l = ; 
Urq = 4(m,+1/m,)v,,+4Z,(m,—1/m,)i,, 
, l . er rie (3.19 
ine = A (m,— 1/m,)v,»+4(m,+ 1/m,)i,, 
= “r 


4. Application to nearly periodic structures 

The resolution of the aperture fields into modes can now be applied t 
the calculation of fields in a structure consisting of a succession of cavities 
not necessarily congruent but having congruent apertures. Every apertur 
except the two terminating ones is the output aperture of one cavity and th 
input aperture of the next. Ifthe input amplitudes for any cavity are know 
the aperture components of the fields can be calculated from (3.5), hene 
the output amplitudes for the preceding cavity from (3.6) and the input 
amplitudes from (3.7). Thus the amplitudes at the first input can b 
expressed as linear combinations of those at the last output; the linea 
equations relating the two sets of amplitudes are then to be solved subject 
to the restrictions imposed by the terminal boundary conditions. If all th: 
cavities are symmetrical the relations between the mode expansions « 
the output of one cavity and at the input of the next can be expressed simply 
in terms of the mode voltages and currents; in fact if functions belonging 
to the first cavity are distinguished by accents, we have from (3.17) and 


3.18) ; 
( vb = Zz Crs Usa? lnd = > di. Yea | 
8s 


8s 


] . 7? Y 1 ‘ , Zz Y 
where Co = | 6,.H#,a8, d,, = +> | 6. 4,48 
h,e, « eh, ry 
S S 


The theory of this section is valid for a structure each of whose cavities 
is of arbitrary shape provided that all the apertures are congruent 4 
stipulated, but its practical value would be small in the absence of further 
simplifying conditions. One such condition is that adjacent cavities should 
be nearly identical, differing from each other by a ‘small’ amount measur: 
able by a small positive number e, so that the whole forms a nearly periodi 
structure and the modes in adjacent cavities are closely related. In this cast 
the coupling constants between different modes in two adjacent cavities 
(the constants c,,, d,, of (4.1) for 7 ~ s in the symmetrical case) will be smal 
quantities of order «, and the methods of perturbation theory will be avai 
able for solution of the equations. If, as is usually the case in practice, tht 
system is operated in such conditions that in zero order approximati0l 
only a single mode and its reciprocal are excited, the modification of the 
amplitudes of the dominant modes due to excitation of the other mode 
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will be of order «®. For calculations of the dominant mode only, accurate to 
frst order in e, the presence of the other modes can be ignored, and the 
astem is reduced effectively to a system of four terminal networks in 
ascade. The v alidity of this approximation is independent of the validity 
{the expansion in transmission modes. The modifications introduced 
by the departures from periodicity are (1) the coefficients in the voltage- 
urent transfer equations of a single cavity ((3.10), (3.19)) have different 
values for different cavities, (2) the output amplitudes or voltages and 
urents of each cavity are not the same as the input amplitudes or voltages 
ndcurrents of the next, but must be obtained from them by further transfer 
equations. In the symmetrical case these reduce by (4.1) to a pure trans- 


former action, ’ — ” ‘ m 
Urb Crp Ura lnb d,, Ura (4.2) 


where, as may be shown without difficulty, c,,d,,—1 is of order e?. 


5, Properties of the transmission modes in symmetrical cavities 

The determination of the transmission modes of a given cavity is a 
boundary value problem of a rather unusual type. Approximate solutions 
can be obtained by standard mathematical techniques for certain simple 
types of cavity, but there is no general method. However, some insight 
into the properties of the modes can be obtained by general methods, par- 
ticularly for symmetrical cavities, and to do this is the aim of this section 
and the next. 

If suitable boundary conditions are imposed on the field at the apertures 
‘,and S, (in addition to the condition nA E = 0 on S’) there will be a set 
f resonance values k,, k,,... of k for which there are fields which satisfy 
these boundary conditions. In particular we consider ‘open-circuit’ 
resonances, with the boundary condition nA H = 0, and ‘closed circuit’ 
resonances, with the condition nA E = 0, on S, and §,. 

Suppose that / is not one of the open-circuit resonance values. Then the 
fields in the cavity, and in particular the aperture electric fields, are com- 
pletely determined by the aperture magnetic fields. The dependence of the 
total field on the aperture magnetic field is linear, and four linear operators 
S... S.», S,,.. and S,, can be defined such that the aperture tangential 
electric fields H, and E, are related to the magnetic fields H, and H, by the 


—" eae G8 &- £84408. (5.1) 


ab 
Since (E*, —H*) is an electromagnetic field in the cavity whenever (E, H) 
is one, the operators S are real (that is, SH is real if H is real). 

If in (3.1) (E,, H,) is taken as the field with aperture magnetic field 
components H, = H,, H, = 0, and (E,, H,) as the field with aperture 








334 E. WILD fLECTE 
components H,, = H}, H, = 0, the resulting equation can be written as [yhere 2 


[ (s 


f.j*d sifticien 
a)’ G0, 


aa aa < 


—_ rx - r ~ 
H,).H3 ds | H,.(§$ er 
S Ss 

which shows that the operator §,,, is adjoint to itself with respect to th: 


scalar product | H,.H* dS. Similarly §,, is adjoint to itself, and §, ; 
J} T#1- 489 Y op J ab | 


adjoint to —§S 

For symmetric cavities the existence of the reflected solution correspond 
ing to any given solution shows that S_,, S,, and $.. Si; thu 
in this case §,,, is adjoint to itself also. 


ba* 


By writing the boundary conditions (2.5) in terms of and the operator 
S it follows that a necessary and sufficient condition that W should b 
the aperture magnetic field of a transmission mode with transmissio 
factor m is that #% should satisfy the equation 





m$,,H#+S,,H4 m(mS,,, 4% + $,, # ). (5.2 
For symmetrical cavities (if the case m 0 is excluded) this reduces to 
S,,% 3(m+1/m)S,,,-%. (5.3 


Equation (5.3) has the same form as the eigenvalue equation for the} The 
resonant frequencies of a mechanical system. The orthogonality relations} transi 
(3.16) can be derived from it directly. If either operator §,,, or §,,, is| more: 


positive definite it also follows that the eigenvalue }(m-+-1/m) must be real 


aa 


and that therefore m must be real or unimodular, as can be seen by multi 
plying (5.3) by #* and integrating over S, 


hence 
| (S,.%).#* dS = i(m+1/m) | (S,,%).#* dS. 
re s but a 
The two integrals are real, and therefore 4(m--1/m) is real unless both thi and 


integrals are zero. 

In the case of a four-terminal network, in which the operators § reduc 
to pure numbers which are necessarily real, it is obvious that (m-+-1/m) ii} This 
always real; (m--1/m) is also real for modes in a cylindrical waveguide. 

The operators, however, are not in general positive definite, and comple 
values of }(m-+-1/m) are therefore possible. An example, for a system wit! 





ess 


a finite number of degrees of freedom, is given by the symmetrical loss 





six terminal transmission network of Fig. 1, in which the quantities X ar 
the reactances of the circuit elements. For this system the operators § cal 
be calculated from the network equations. The equation for the eigenvaltt 16. C 
A }(m-+1/m) is found to be In 
A?(XF-+ a?2.X3)— 2Aa(1—a)X2+ {(1—a)?X2—X'! 


2) | 
if 0, in ¢ 
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tenas [Pybere x = X,/(2X,+X,). This equation obviously has complex roots for 
sificiently small values of a (small values of X,) if X} > X} > 0. In fact 
the roots are complex if 

2a)X3 > X? > 0. 
t to th 
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id S.,, ! = i 4 1 
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ces 
Fic. 1 
for t The relation between the aperture electric and magnetic fields for a 
elations| transmission mode in a symmetrical cavity can be put into a form which is 
r §,, is) more convenient for later developments. We have for such a mode 
he ~ ~ 
iE, Lmd mS a H+ S., HK , 
ym 
, o, | Z 
ence uc Sac H T Sa H ’ 
m 
ut also iE, id mS,,, a A 
oth . 
nd therefore from the relations oo au. >. a. 
> red id }(m—1/m)S8,,,%. (5.4) 
) Tl ° . . ° ° 
this gives an alternative formula for the normalizing integral 
rude 
omplex} eh, 6... HAS 
mM Wi . 
lossles , 
, (a | (S,, 4%). 4%, dS8. (5.5) 
sad 2 mW : . , 
's S S 


6. Operation at and near resonance 


| Ina symm trical cavity, since the reflected solution corresponding to 


? 


open-circuit or closed-circuit resonant mode is another such mode, 
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these modes are, or in case of degeneracy a linear basis for them can be 
taken to be, all even or odd (that is, to have the non-vanishing aperture 
tangential components at the input and output apertures equal or equal 
and opposite). This being the case, such a resonant mode is also a transmis. 
sion mode with transmission factor +1 for an even mode, —1 for an odd 
mode. Conversely, if there is a transmission mode with factor m = +] 
which is neither an open-circuit nor a closed-circuit resonant mode the sum 
and the difference of this mode and its corresponding reflected solution are 
open-circuit and closed-circuit resonant modes, even for m +-1, odd for 
m —1. Thus transmission factors +1 occur at and only at resonant 
frequencies, and we wish to study the behaviour of these modes and of the 
modes which approach them in the neighbourhood of resonance. 

Consider then a simple (non-degenerate) closed circuit even resonant 
mode. Two cases arise. 

1. There is also an open-circuit even-resonant mode of the same fre- 
quency, and the aperture fields of the two modes are not orthogonal, 


| E.HdS +0. Then any linear combination of the two modes is a trans- 
S 


mission mode with factor +1; such a linear combination and the corre- 
sponding reflected solution form a pair of reciprocal modes to which the 
theory of section 3 can be applied directly. 

2. There is no such open circuit mode. This case requires further con- 
sideration. 

In case 2, let H, be the tangential magnetic field in the input and output 
apertures for (arbitrarily chosen) unit amplitude of the resonant mode. 
Since there is no even open-circuit resonant mode there is a (unique) 
electromagnetic field in the cavity which has boundary values }% and 

1H, for the tangential magnetic field at the input and output apertures 
respectively, and by symmetry the tangential electric field for this mode 
will be the same, say &, at the input and output apertures. This field and 
the resonant mode obviously form a pair of associated solutions of the type 
described in equations (2.6). Since the resonant mode is a transmission 
mode with factor 1, % is orthogonal on S to the aperture tangential com- 
ponent &, (r = 1, 2,...) of the electrical field of every other transmission 
mode; and by applying (3.1) to the field of such a mode and the field with 
boundary values (&, }%) and (6), —3%) one can see that & is orthogonal 
to every aperture field #. Thus the sets (&, &,.-.) and (%, #,..-) are 
bi-orthogonal, and we assume them to be complete. (The set (6, &,-.-) i 
obviously not complete since %, is orthogonal to every member of it.) 


We now put [ &,.%,d8 = 1/Y, (6.1) 
s 
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idefine the voltage vy and current 7, of the Oth mode in a given electro- 
netic field (E, H) at the aperture S by 
Y | &.Ads, (6.2) 
S re 
» that the expansions (3.17) hold with hy = 1 and e, = 1/Y. 


The transfer equations from output to input now become 


Vo, “0,b> lo,a Y vo, T lop: (6.3) 

Equations (6.3) are the limiting forms of equations (3.19), with r 0, 

tends to 1 and Z, tends to 0 in such a way that }(m,—1/m,)/Z, tends 

| to. It will now be shown, with reasonable assumptions, that this limiting 
wess is significant in that, for small variations of the shape and size of 
ecavity which involve a departure from resonance, there is a pair of 


iprocal transmission modes whose voltage and current relations tend 





ythly to those of (6.2) and (6.3) as resonance is approached. It is 
sumed for this purpose that there are no open-circuit resonant modes 
/ofeither parity at the operating frequency. 

The aperture field W and the eigenvalue A }(m-+-1/m) for a transmis- 

n mode are determined by equation (5.3) 

S,,% =)S,,, 
the ' indicating functiuns corresponding to the distorted cavity) and the 
perators S),, and §,,, may be assumed to vary smoothly with the dimen- 

sions of the cavity since the resonant mode involved is a closed circuit, 
| not an open circuit. For the undistorted cavity the equation has a non- 
legenerate solution # = #, with A = 1, and it may be assumed that for 
uall changes in dimensions there will be a solution “4, with eigenvalue 2’, 
which by suitable choice of the arbitrary amplitude factor can be taken to 
| tend smoothly to #9, A’ tending smoothly to 1, as the changes in dimensions 
tend to zero. This solution gives rise to two reciprocal modes with the two 
solutions mj), and 1/mj), of the equation $(m-+-1/m) \’ as transmission 


> 


ctors. We take for these modes hj = 1 so that e, = Z, and, from (5.5), 


: ; | , > | am 2 —re 
Z ( =: (Mo—1/mp) (S,,% 0). 4%, 48, 


S 

| While Irom (5.4) the aperture electric field & is given by 
| es ll 0) Sia %% 

. Lye mo) ba 7 9 
a= t 
| In the undistorted cay ity there can exist an electromagnetic field with 
soundary values (&, %)) at the input aperture and (é, 0) at the output 
tf onartiarn ° . P i 
| ‘perture, this being the sum of the resonant mode, with boundary values 
Z 
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(0,3.4%,), (0, 4% ), and the field with boundary values (Eo, Ho), (6, —4H, 
Thus in the undistorted cavity S,,% = Spa %o = 16, and S$), H;, tend 
smoothly to ié, as the distortion tends to zero. 


Thus the functions a , 
$(mo+1/mo), 


, , r7t l 9 i , we! gz Y , , 4 ad 
$(Mo—1/my)Zy = ; (A’2—1) | (Sjg%o)-4%o dS, 5 (29 — 1/m)/Zy, 


S 

KH y/hy and &)/e, tend smoothly to the limits 1, 0, ¥, %, and Y & respectively 
the coefficients of the transfer equations and the voltage and current of 
fixed aperture field tend smoothly to their values for the undistorted cavity 

In general mj, itself and $(m,—1/m 5) will not vary smoothly since ; 
}(mo+1/m) = 1+ with € small, my ~ 1+,/(2e). If. is positive my is rea 
if « is negative m, is complex and unimodular. The above results apply t 
variations of the operating frequency as well as to distortions of the cavity 
so that in general a resonance of this type separates a pass band from ai 
attenuation band of the mode concerned. 








The methods of this section can be applied to an odd closed-circuit 





resonance, and the results hold with small modifications in detail. For th 

























discussion of open-circuit resonances a development parallel to that oj 
section 5 with the roles of the electric and magnetic fields interchanged i 
necessary. Smooth variation of the voltage and current relations is ther 
obtained by taking e, = 1 so that hj tends to 0. 
The phenomenon of an open-circuit and a closed-circuit resonant mock| 
of the same parity at the same frequency (case 1 above) occurs in a structure 
of n identical symmetrical cavities (considered as a single cavity) whe 
sach cavity has a transmission mode of factor m such that m” +1. Sue! 
a value for the individual cavity is of no special significance and this suggest: 
that case 1 in general is an accidental degeneracy of no particular im 
portance. The transmission modes in the neighbourhood of such an eve 
resonance can be shown to be determined by an eigenvalue equation 
S,¢ =i HK, $,.# = i= é, 
m+ 1 ‘i m--] 
where S, = S,,,—$ 


give the aperture magnetic fields as functions of the electric fields. Ti 


sa Sp, and §, is similarly defined in terms of operators whit 


eigenvalue m varies smoothly in the neighbourhood of m 1 and for an 
one-parameter family of variations of $, and §, there is a pair of transmis 
sion modes which tend to definite linear combinations of the open-ciretl! 
and closed-circuit resonant modes. This is the behaviour which would ! 
expected of an n-cavity structure of the type mentioned when the operatiti 
frequency is varied. 
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* aI 7, Application to a practical structure 

sia As an example of the application of the method developed above, con- 
sider a structure of the type used in the Berkeley proton accelerator (2). 
This is developed from the periodic structure formed by the repetition of 

» dements sketched (in somewhat idealized form) in Fig. 2. The protons 
A 

tive | es 

‘ent 01 | Ff 

1 CaVl | | 

Sit ce f 
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For t | le § | 
. that asians emai. 


anged N A 
s is the Fic. 2 
int mode} travel along the axis of symmetry OO’ and are accelerated by the electric 


structut} field in each gap BB’. The length L, which is determined by the relation 
ty) whet} LA = vc (A the operating wavelength, v the velocity of the protons, c the 


1. Sucl| velocity of light) increases slowly along the structure to allow for the 
s suggests} acceleration of the protons; each cavity is designed to be a closed-circuit 


cular it | resonant cavity by making the dimensions g and d vary with L, the outside 
, an evel} diameter D being constant. (The resonant mode is one in which the only 
tion | on-vanishing field components are the cylindrical polar components 


i, E,, and Hy.) In the Berkeley accelerator this variation is achieved in 


the early stages by increasing g/L, keeping d constant, and later by in- 


teasing d keeping g/L constant. For simplicity we consider only the former 





ail 
oo T | type of variation. (If d is varied it is necessary to take the dividing plane 
3 a between cavities as the central plane AA’ of Fig. 2 in order to have congruent 
transmis pertures.) If each resonator is accurately tuned, the input-output voltage 
i ind current relations are given by (6.3), and if the accelerator is terminated 
esa by a short circuit the voltage is zero at all the apertures. The variation of 
operatit irrent (and hence of the electromagnetic field strength) will be determined 


j °Y the transformer action at the junction of adjoining cavities ((4.1) and 


t.<)), and it is this variation that we calculate. 
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The resonant field in the cavity can be calculated approximately by 


assuming that the longitudinal component of the electric field on the 


cylinder BB’ which forms the continuation of the inner conductor is 
constant. The field in the space between the coaxial cylinders can be 
obtained as the sum of a series of coaxial waveguide modes, and the aperture 
field # is thus determined. In a similar approximation the corresponding 
‘odd-magnetic’ field of section 6, which determines the aperture field ¢, 
can be calculated assuming zero longitudinal electric field on BB’. The 
integrals (&.# dS and | &.#' dS, which determine the ratio of currents 
Ss Ss 

in adjacent cavities, can be obtained conveniently from these series in 
virtue of the orthogonality relations of the waveguide modes. 

It is convenient to define unit amplitude of the resonant mode in the 
cavity as that which corresponds to unit value of the mean longitudinal 
electric field on BB’.. With this definition %, is found to have the value 


ik . 7K, oe d 2)Jg(k,, D 2) sinh(m,, g 2) “a - 
Zl mn Wolken d/2)?—{J,(«,, D/2)?? sinh(m,, L/2) Aen?) 
l sin(kg/2) 1 


~ klog(D/d) sin(kL/2) r|’ 
where «,, K,... are the roots of the equation 
Jo(«d/2)¥o(«D/2)—Yo(«d/2)Jy(«D/2) = 0, 


m \/(«2 —k?), Jy and Y, are the Bessel functions of first and second kind of 


n 


order zero, and 


K(k, 1) = Yj(x, D/2)J,(«,, 7) —Jp(x,, D/2)¥;(«,, 7). 


i 


é, is 
a ~_ Ke /2 2 ; 
ol Jo(n d/: 2)Jo(%n D/2) , sinh(m, 9/2) coth(m,, L/2) f(x 
21 m, {Jy(«,, d/2)\?—{Jy(«,, D/2)\? sinh(m,, L/2) 


; __! sin(tg/2) cot(kl 2)! 
k log(D/d) sin(kL/2) r 


The ratio i’/i | E,H,d8 || €,#4dS of currents in adjacent sections 
Ss Ss 

in one of which L is replaced by L’ = L+8L and g by g’ = g+8g, can bt 

found to first order in 5L and 8g by evaluating the integrals and differ 


entiating with respect to L’ and g’. The result can be expressed in the form 


dlogi = (1+ F,)6L/L—(1+- F,)dg/g 
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ynere 
> G,{(m,, L/2)coth(m, L/2)—1}4Go{(kL/2)cot(kL/2)—1} 
F,. = - ¥ scl acmnmaccsesindiacaiscds 
> G, T Go 
n=1 
> G,{(m,, g/2)coth(m, g/2)—1}-4 Gf (kg/2)cot(kg/2)—1} 
F ) 1 inne 
> G,, T Go 
n=1 
{J,(«,, D/2)\? sinh?(m,, g/2 
G. ’ ol, ) " — nJ ‘ ) coth(m,, L/2), 
m? {Jo(«,, d/2)}*—{Jy(«, D/2)}* sinh?(m, L/2) 
l sin?(kq/2) 


2 log(D/d) sin*{ kL )2) oer "/”)- 
Ff, and F, are small compared with unity for the dimensions considered. 
Ifthey are neglected the equation for i can be integrated and shows that 
) Listhe same for all cavities. ig/L is proportional to the local mean value 
E of the longitudinal electric field along the axis of the system, and there- 
fore to the velocity increment per cavity. The design of the Berkeley 
accelerator was based on the assumption of constancy of this quantity. 
Values of F,; and F, were computed for the range of values 
0-09 <= L/A < O18, 
with the values D/A = 0-66, d/D = 0-124, used in the Berkeley accelerator. 
lhe values of g for resonance were obtained from the empirical formula 
given by Alvarez (2), the range of g/L being from about 0-2 to 0-3. In this 
range it was found that F, was negligible (< 0-005) and that F;, could be 
ccurately represented by the formula F;, 0-64(L/D)?. With this 
formula the equation for i shows that EF is proportional to ¢—%3%4/D", 
giving a variation of about 2 per cent over the range of values considered. 
\ variation of this amount would be completely masked by tuning errors, 
but the form of dependence on L suggests that the effect might be important 


t the high velocity end of the accelerator. 
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ON A CLASS OF DIFFERENTIAL EQUATION 
GOVERNING NON-LINEAR VIBRATIONS 


By A. W. GILLIES (Northampton Polytechnic, London) 


[Received 4 October 1956] 


SUMMARY 
A differential equation of the form 
f(D)x+g(D)y = Fit) 

is considered, with y a function of x, expansible in a series in powers of x; F(t) isa 
simple periodic function and D is the operator d/dt. The approximate linear equatior 
when x and y are small is assumed to have one nearly sinusoidal normal mode, its 
other modes being exponentially damped with time constants less than a constant 
which is O(1). 

The method previously applied to particular second- and third-order equations of 
this form is applied to determine non-resonant and resonant periodic solutions of this 
equation. The solutions considered are of the form 26 cos(wt+¢) plus terms of higher 
order, where } and ¢ are either constants satisfying an amplitude equation or approxi- 
mate to slowly varying functions of the time, satisfying an autonomous system of 
variational equations. These are of the same form as in the particular cases. Thus, 
although the equation is of arbitrarily high order, it behaves, as far as the solutions 


considered in this paper are concerned, like one of the second order. Although th 





solution of an equation of order n depends in general on n arbitrary constants, it 
this case, in consequence of the form assumed for the linearized equation, only the tw 
parameters 6 and ¢d are of importance, apart from a transient stage, the duration of 
which is of the order of the longest time constant of the linearized equation. 

The variational equations are of four main types; in two of these the final state is 
a periodic solution or a combination oscillation determined by a singular point or 
limit cycle of the variational equations, in the other two the oscillation amplitud 
may increase until the variational equations cease to be valid, and the final stat 
cannot be determined by the methods of this paper. 


1. Introduction 
In a large class of oscillation phenomena one is concerned with a system 
which is substantially linear with the exception of one non-linear connexion, 
which leads to a differential equation of the form 
f(D)a+q(D)y = Fit), (I 
f(D) and g(D) being polynomials in the operator D = d/dt, with a non-lineat 
functional relation between the dependent variables 2 and y. If this relation 
is analytic then for sufficiently small values of the variables a power series 
expansion of the form 
9 ae 9 
Y = €, 4+, par?+C, 723+... \4 
in which c, may be assumed positive (otherwise y could be replaced by —Y) 


(Quart. Journ. Mech, and Applied Math., Vol. X, Pt. 3 (1957)] 














the ¢, i 
charact 
verge 1 
In p 
ften | 
an eXf 
will tl 
ippror 
Ine 
to an 
terms 
the fir 
Wit 
forein: 


The 
pheno 
has al 


It i 
can b 
The | 


with 
this 1 
in an 

By 


perio 


with 

It 
time 
Le. a 
plan 
from 

( 
the | 




























ON A CLASS OF DIFFERENTIAL EQUATION 343 


thec, in general are assumed to be O(1) and yp is a small positive parameter 
haracterizing the magnitude of the non-linear terms. The series will con- 

erge in some interval p|x k,. 

In practice the analytic form of the relationship between 2 and y will 
ften be unknown and will be replaced by a polynomial approximation to 
n experimentally determined relationship. The interval of convergence 
vill then be replaced by the interval of validity of this polynomial 
pproximation. 

Ineither case it will be assumed, not only that the variables are restricted 

in interval in which (2) is valid, but that furthermore the non-linear 





“(t) isa } terms are small, so that in the equations of the first approximation only 


the first three terms of (2) are used. 





With x small enough for non-linear terms to be neglected, and with the 


forcing term F(t) also omitted, we have the approximate linear equation 








| {(D)+ce,g(D)}x = 0. (3) 

t high The possibility of free oscillation, and of resonance and synchronization 

PPT phenomena with a periodic forcing term, will arise when this linear equation 
Thus is an oscillatory normal mode which is at most only slightly damped. 
tior It is therefore assumed that the polynomial 

we G(z) = f(z) +e, 9(2) 

in be factorized in the form 
, 
} G(z) {(z—e)?+ O7h(z). 
statels | The linear equation (3) then has a normal mode of the form 
Ae cos(Qt-+-d) 

_ | with A and ¢ arbitrary constants, and if € is small in comparison with Q 
| this will approximate to a simple periodic function, the fractional change 
| namplitude in one period being small in comparison with unity. 

ois | By a change of time scale the equation may be normalized to give the 

cat period 27 for this normal mode, when we shall have 

G(z) {(z—e)?+ 1 }h(z) 
| with « small in comparison with unity. 
linea! It will be further assumed that all other modes of (3) are damped with 


Jation | “Me constant less than 7’, where 7’ is a constant not greater than O(1); 
e. all the zeros of A(z) are to the left of a line rez 1/T on the complex 
plane, fA(2 will then have a minimum value, say m > 0, as w varies 


irom to f 





(ww) will then have a minimum value greater than m/2 if w? is outside 


ne interval w* < $ around the resonant frequency value 1. 
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it is non-singular. 


2B cos wt. 


F(t) 


{ f(D)+ 


nn 


ce, g(D)}x+g(D) ¥ c, peta" 


It is also assumed that the degree of g(z) is not greater than that of G(2). 
i.e. that the order of the original equation (1) is the same as that of the 


linear equation (3). In other words, regarded as a perturbation problem, 


Finally, the forcing term is taken as a simple periodic function 


Thus, on eliminating y we have an equation of the form 


2Becoswt 


{(D—e)?+1}h(D)x+9(D) ¥ c, pw" a" 2B cos wt, (4) 


n” — 


in which 
i. 
IT. 


e and pw are small parameters with p > 0; 


min |h(iw) m>0O for —wo<w< 


from which follows 


min {(iw—e)?+1}h(iw)| > dm 


the zeros of A(z) are all to the left of rez 
not greater than O(1); 


the degree of g(z) is not greater than that of 


f(2)+e,9(2z) = {(z—-6)? + RA(z). 


in when required. 


that 


of variational equations [(37) and (38) below]. These equations, as We 


+; 


for 0 < w? 


9 


} and 3 < w?; 


—1/T with T' a constant 


It is the purpose of the present paper to consider the periodic solutions 
of (4) subject to conditions I to IV with particular reference to the resonance 
case when w is near to 1, i.e. w—1 O(e), B is O(e), and the small para- 


yu”. These latter assumptions will be brought 


In previous papers, particular second order (1) and third order (2 
equations of this general type have, been discussed. The same procedur 
is applied to (4) and results of the same form are obtained. Although th 
equation is of arbitrarily high order, and the solution is dependent on 
corresponding number of arbitrary constants, it appears that, owing t 
the form assumed for the linearized equation, after a transient interval i 
as the solutions considered in the paper are concerned 
equation of the second order. These solutions are of thi 
¢) plus terms of higher order, where ) and ¢ are eithe 
satisfy an amplitude equation [(17) and (18) below] & 


slowly varying functions of the time that satisfy an autonomous syste! 
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ON A CLASS OF DIFFERENTIAL EQUATION 345 
sthe resonance curves, are of the same forms as were obtained in the 
articular cases. 

?, The non-resonant solution 
We seek to determine a periodic solution of (4) in powers of B. Set 
x YD B+ al?) B24 73) B34... (5) 
Substitute in (4) and determine the 2” in succession as periodic functions 
ft by equating coefficients of like powers of B. 


f(z) +e; 9(2) 


Writing C(z) 
g(z) 
nd # C(riw), 
the first few a” are yD e+e, 
th 
at l all it | 


f(tw)- Cy g(tw) 


vc), = conjugate of a{) ; 
x re x) 1 ones (6) 
th 9 
9 Co px,’ 
a — = P 
2 t 
2c, pra) 
mS = > 
=a - 
So 
) , — 
p'?), conj. of x; 
(3) (3) _| »(3) | (3) 
A 3 44 L113, 
vith « 
; is 3 
2Cok 2 rs- Cy ph y(1)8 c2 pta\)) 
ye = “ - — Ce ~ ; (7) 
$3 \So G3 
¢ yal yD) r\?)) -3¢ L y(1)? (1) 
(3 Ve PrA~) 1“2 it aded 1 
l ¥ 
$1 
1¢2 Ie2 2,.(1)2 (1) 
2, “2 «o, Kh *-1 
Ole, a e 
So . $1 


ind x's, 2° the conjugates of these. The general form is 


*summation being over alternate integers from —n to +n. The term 
contains the exponential factor e"! and combines with its conjugate 
to lorm a real term of rth harmonic frequency. 

With these values for the x” the series (5) gives a formal solution of the 


itte . } . ° 
lilferential equation which is of period 27. 
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It may be noted that 2” contains »”-! as a factor, so that the series 
may be regarded as a series in powers of p. 


3. Convergence of the non-resonant solution 

The absolute and uniform convergence with respect to ¢ of the foregoing 
solution may be established by the aid of the equation 

E = Nyn—cp, n?—cp? 7 —... 

in which c is a positive constant greater than |c,,| (n = 2, 3,...) and N isa 
positive constant to be chosen later.+ 

This is equivalent to 

cy 
] — Ly ” 


which defines 7 as a two-valued function of €. The branch which vanishes 


& (N- c)n— 


for £ = 0 may be expanded in powers of € in a circle extending to the 
branch point nearer to the origin. This is at € = R, » = S, where 


R=- [N+ 2c—2(Ne+c?)4], 
My 


1 
Pen |e ( Poe. (3 
By N+e} | 


If this series for 7 in powers of € is obtained by reversion of the former 


series it will be seen to be a majorante of (5) provided 


BS Py, (9 
__—2B = é os RK (10 
Flin) Fe, giw)| ~ NN 
C|>N (r= 0, +1, +2....). (11 
In fact, if the series is written 
7 ADE + nM? + ME +... (19 
and (9), (10), (11) are satisfied, we have 
)| Br ; > ia < Fs. 


The ¢ sries (5) will then converge absolutely and uniformly with respect 
to t, provided we : 
} R f(iw)+e, g(tw) 13 

N 


for any » < p,; and it will then provide a periodic solution of the differen: 


B< B, 


tial equation (4). 


+ In particular cases a wider interval of convergence may be obtained using 
g Nn — |€o|419?— |C3|41 7° — «++ 


but this is less convenient for the general discussion as it cannot be expressed in finite 10m 








Ife 


mined 


provic 
suffici 

If, | 
is unr 
satisf 
{(z) d 


series 


4, Tk 

Wh 
be ust 
small 
If thi 
Y inc 
will b 
(.; a 

We 


2 COs 


with 
wher 
Thus 
Su 
tainit 
giver 
omis: 
the h 
Th 


form 
en tot 
in e 
equa 
term 


f( dD 


or 


€ series 


reg iT g 


LN is 


anishes 
r to the 


e 


- forme! 


ifterel 























ON A CLASS OF DIFFERENTIAL EQUATION 347 
Ife is independent of », then B, may be chosen arbitrarily, and py, deter- 
nined by (13). We then have convergence in any specified interval of B, 
yovided is sufficiently small, ie. provided the non-linear terms are 
sficiently small. 

If, however, € is small this may require that p is smaller than O(e?) if w 
sunrestricted. This arises from the fact that N must be taken O(e) to 
tisfy (10) if any rw is near to 1. [Note that condition IV ensures that 
‘:) does not tend to zero at infinity.] That is, if resonance occurs the 


eries solution (5) will only have an interval of convergence of order «?/y. 


4, The resonant solution 

When an rw is near to | the solution of the previous sections ceases to 
e useful when € is small. The failure arises through the occurrence of a 
small value, of O(e), for one of the ¢, divisors in the determination of the 2”. 
Ifthis particular divisor can be avoided, it may again be possible to choose 
Y independently of « and so secure a useful interval of convergence. This 
ill be done for the case of fundamental resonance when « is near to 1 and 
(,, are O(« 


We suppose that a periodic solution exists whose fundamental is 


%eos(wt+d) and set 
x DDH + aH? +- aH3 +... (14) 
h a e+ a), 
vhere yi ‘ +p) a! S ¢ i(ot+¢) 
Thus 2h 2b cos(wt+¢) is the complete fundamental. 


Substitute (14) into (4) and choose the 2” so as to annul all terms con- 
taining 6” other than those of fundamental frequency. The 2” are then 


given in terms of a{? and a2, by the same expressions as before with the 


mission of wv and a) when nv is odd, and of terms arising from these in 
the higher orders. 
The terms which were previously cancelled by the 2) now remain to 


lm an equation containing only terms with e‘ and their conjugates with 
his equation will be satisfied if the equation formed by the terms 


‘ only is satisfied; then 2® will be given by (6) and the following 
jwations, and 2 will now contain only x as given by (7) and 2°). The 
terms in which remain are 
2 Dp2 \ 
f(D) +c, g(D) 2 b—g(D i = Th Be) ? xl) 534 a Bei (15) 
Co So 


C1 g9(iw) tb g(iw){ =? + = — Be, )u2b? Lo. = Be id. (16) 
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Equation (16), equivalent to two real equations, is the amplitud 
equation, and if values of b and ¢ exist which satisfy this equation and 
for which the series (14) converges, then each such pair will determine , 
periodic solution of the differential equation (4). 

By this procedure ¢, , do not occur as divisors, so N has to satisfy (10) onh 


for r= 0, +2, +3,... and not for r= 1. If w is near to 1 and nox, 
(r = 0, +2, +3...) falls within 4 < r?w? < 3, then by II it is possible t 
choose N independently of «; (12) is now a majorante of (14) provided 
(9), (11) with r + 1, are satisfied, and instead of (10) 

b< R/2N 


with R still given by (8). 

We can now choose 6,, then p,, so that b, < 4R/N and we have that 
(14) converges absolutely and uniformly with respect to ¢ for b < b, and 
re 

Furthermore, the sum of the moduli of the terms containing b” in th: 
amplitude equation is less than NV7"€" so that (12) multiplied by JX is 
majorante of the series in the amplitude equation (16). This therefor 
converges in the same interval. 

The determination of periodic solutions of period 27/w and amplitud 
less than b, is now reduced to the determination of solutions of the ampli: 
tude equation. 

It may be noted that 2” in (14) contains the factor «”~!, so that (14) may 
be regarded as a series in powers of ». Likewise the amplitude equatia 
has the general form 


Be id <A 
‘ N A ane +1 
q( vw ) n=0 


’ 


where A, = ¢, and A, (n > 1) is a function of the coefficients ¢,...., ¢, 

and of the ¢_, for r= 0, +2, +3...., +2n. Thus this equation also ma 
be regarded as a series in powers of x. In fact this viewpoint is appropriat 
in considering the orders of magnitude of the terms, once we have fixe 


the interval of convergence b < b,. 


5. The solution of the amplitude equation 

We know already that if B is sufficiently small and €, = C(iw) is not zer 
the amplitude equation has one solution, given by the non-resonant so! 
tion, for which 6 > 0 with B. The non-resonant solution could in fact! 
recovered by reversion of (16). 

To consider further the solution of the amplitude equation, we mls 
make some assumptions regarding the magnitude of its terms. We thet 
fore introduce the following: 
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9 


V. the small parameters are related by |e| = p?; 
VI. the amplitude of the forcing term is O(e) with B = |\e\E; 
VII. the frequency w differs from 1 by O(e) with w—1 = |e\c. 


With these three assumptions we have 


f(t) +e, 9(%) 2ih(i)e+ O(c?) = eAe’*+ O(e?) 
with A 2\h(t)\, x = argh(i)—47; 
f’(i) +e, 9'(%) 2ih(i)+ O(e) = —Ae’*+ Ole). 
So f(iw)+c, g(tw) Ae’™(§—ia) \e| + Ole?) 
ac §=5=+41. 
€ 


he factor occurring in the coefficient of b? in (16) is 

i fee tes 2c , os ema 
: =e = Be —_ =. _ 3¢,+ O(e) = A’+ iA" + Ole), 

=v) 
say, and g(iw) g(t) T Ole). 
[he amplitude equation now becomes 
Ae'*(6—ia)b—g(t)(A’ +7A" )b3 + O(n?) Ee-'?, 
Omitting the O(u?) we have the equation of the first approximation, which 
may be rearranged as 
g(r)(A —T iA ) p3 ET ib x), 

Ae’ A 

g(i)(A’ +iX") 


Ae 


{O 1a0)b 


A eid 
and separate real and imaginary parts. Then 


d5b— A cosa 68 - contd +x), (17) 


: , E . 

ob+AsinaAb  Sin(e +), (18) 

which gives the real equations of the first approximation for b and ¢. 
When the linear equation has an increasing oscillatory mode e is positive 


and § |. In this case (17), with E = 0, gives for the free oscillation 
, | 
b 0 or 62 > 
A cosa 


Whust (18) gives for the second value of b 
o A sind b? —tand. 


‘or a free oscillation to exist, cosA must be positive and not small of O(y) 
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for the value given for 6 to fall within the interval of convergence b < }, 
For this to be the case we require that 
VIII. —$7+A <A<}n—A withe>0. 

When this condition is satisfied, let the free oscillation amplitude be a: 


(17), (18) may then be rearranged as 


“(! ~ =) E cos(f-+«) 


a a2 Aa 
b E sin(d+ x) 
“(e+ a —_— 
where y = tana. 


If we write E/(Aa) F,d+a= ¢’, b/a = b’, these equations become 

b’(1—b’?) F cos ¢’, (19 

b'(o-+vb’?) F sind’. (20 
These are of the same standard form obtained in the particular cases (1), (2), 
The solution of these equations has been considered previously and typical 
resonance curves (graphs of 6? against o for fixed F’), are given in (1), 
Figs. 1 and 2, for vy = 1 and v = 2, respectively. For v = 0 they are the 
symmetrical van der Pol resonance curves. The curves for v < 0 are 
obtained by reflecting those for v > 0 in the b? axis. 


These resonance curves are therefore typical of a non-linear oscillatory | 


system having a nearly sinusoidal free oscillation. 

We may use implicit function theory to show that the complete equations 
for b, ¢ have solutions when |e pu? is sufficiently small, differing at most 
by O(e) from the solutions of the equations of the first approximation. 


6. Differential equations for the amplitude and phase 
Denote the formal solution (14) by 
K(b, d, t) = 2™b+ xb? + 23+... 
2b cos(wt+¢)+ O(n). (21 
When K(b,¢,t) is substituted for x in (4) the residual terms form the 
amplitude equation. 
Write the differential equation (4) in the form 
(D?+1)h(D)x—(2eD—e*)h(D)x+9(D) ¥ ¢, p"a"—2Beoswt = 0 (2 


n 
or (D?+-1)h(D)x+-pS(x)—2B coswt = 0, 
in which pS(x) denotes the complete group of small terms, and is actually 
a function of x and its differential coefficients, i.e. 


peS(2) — (2u?D—p')h(D)x+9(D) 3 jw". 
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The residual terms on the left when K is substituted for x are 


4ce2 2c? re ' ‘ i 
- < Jn 3c) 2U+. ef _ Beivt + con}: 
So So 


. 9 4H (b, ¢, t) 


Let these terms be denoted by 2 7 . 
tf 


ee : 
1} f(ww) C, g(iw)}b q(tw) 
| 


} 


wt @. 





in which @ 


H(b. d, €) [eres f(tw)- 


2tw 


C, g(iw))b+...} — Be’]+-conj. 


12 

f - [e'"} Ae'™(d 
») 
at 


ia)b—g(i)(A’ + 1A" )b3}— Ee’+- O(u?)]+-conj. 


fh. d are variables along with t we have 


aa I |h{— | K uwS(K)—2B cos wt = 


9 
c 


’ , | » ¢ C r c . 
uS(K) (2p>— wih) K t <) a 


n 


with e. g*-1K*. 


The differential equation (22) is now replaced by the system of two 


equations 


(D2-+-1)u 2Becoswt 0, (24) 


(25) 


pS(2) 


w= MD)e; 


the formal solution of which is 
v kK (b, ¢,t), 
u — h(D)K(b,4,t) 
K,(0, 4, t), 
say, with b, dé constants satisfying the amplitude equation, i.e. for which 
H and dH dt both vanish. Now take b, ¢ as new variables by using 


u kK, (5, 4, t), (26) 
du ok, H. (27) 
dt ct 
rm c 4 db 9 ¢ 
These require that Ay t | OK, dp | H 0. (28) 
; cb dt Ch dt 
Substituting into (24), using (23) and (25), we obtain 
c ( K, db C jc K dd cH 
(eee ES [—_> _ A 1 pt hie) — Sik 4 — = @. 29 
Aer - al oy a Hy (x) (K)} Ot (29) 
Solving (28) and (29) for db/dt and dd/dt, 
db 1(_oK 0 (ek | 
a ( H) 30 
dt A,\ Ch Ch\ ct j we 
dd 1 | oK 0 (eK | 
L—+.8 +}, 31 
dt A,\ eb ail ot } (31) 
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where L = p{S(x)—S(K)}4+ ” 
C 
OK, 0 (ekK ok, e(eK 
' Bh: on See ee fe « '_H). 
an 1 3b ad\ at H) ad | at 


Instead of the original equation we now have the system of three equations 
(30), (31) together with 
h(D)x K,(b, 4, t) (32) 


in the three unknowns b, ¢, x. These are satisfied if b, d have constant 
values satisfying the amplitude equation, and x = K(b,¢,t), i.e. by the 
periodic solutions of the original equation. 

The right-hand members of (30), (31) may be expanded in powers of u, 
Evaluating only the terms of lowest index we obtain 


® _# £S(a) —S(K)}cos(0+ x) +O(p)}4 
dt A 
+p3)3b- A cosaA 6? — = cos(¢- | +O(p3), (33 
ye = ; {S(x)— S(K)}{sin(@+ «)+ O(n)}— 
tf Z 


— {ob AsinA b? — : sin(d-+- 0} + O(u3). (34) 


Under VIII, when a free oscillation is possible, these become with the 
notation of (19), (20) 


db’ BM to ae 
ome ,S(x =f ‘Seos t 4 gee: Soden 
hs (x) —S(K)}{cos(wt+¢')+ O(p)} 
+-p*{b'(1—b’?)— F cos ¢'}+ O(p'), 
yc i < 7 {S(x)—S(K)\{sin(wt+4')+ O(n)} 


—p*{b'(o+-vb’?)— F sin ¢’}+ O(u'), 


which are of the form discussed in (2). 


7. Qualitative character of the solutions 

If the differential equation (4) is of order n*, a solution will be determined 
if values are given for x, Dz,..., D”*-1y at t = 0. Such a solution may be 
represented by a trajectory [in a space Y in which 2, Dz,..., D”’~x are 
cartesian coordinates; the trajectory starting from a point P, whose 00 
ordinates are the initial values. 





Th 


hz) | 

A’ 
b,o1 

Th 
inter 
y+ | 
of F 
unif 

W 
of F 
the « 
R' a 
to (b 

T 
at P’ 
for 
appr 
tion 

Sit 
than 
remé 
O(1) 


Le, 


whe} 
coeff 
resp 

Le 
to tl 


ations 





ith the | 


rmined 


be 


may 





the cylinder b < b, for all t > 0. Any point P in Z transforms to P’ in 
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The space ¥ is transformed into #’ by the non-singular transformation 
H i x 
Dr 3y D- 3y 
u h(D)x 
du/dt = h(D)Dzx, 
being of degree n 2: T will transform into I’. 
A further transformation to “” is made by (26), (27), polar coordinates 
S replacing w and du/dt, the other coordinates being unchanged. 


The latter transformation is only meaningful provided b < b,. To the 


nterior of the circle b < b, in the b, d plane of Y” corresponds a region 
lu/dt)? Ab,—O(y) approximating to a circle in the u, du/dt plane 


/’. and this independently of t since the series in (26), (27) converge 
niformly with respect to ¢ 

We can therefore certainly define a sphere # with centre at the origin 
f¥Y which transforms into # in’ and &” in. FY” such that Z” is within 


and P” in 42”. The determinant of the transformation from (u,du/dt) 
)is A, which vanishes only when 6b = 0 and ¢ is indeterminate. 

lo a trajectory [ starting at P, in Z correspond trajectories I’ starting 
t P)in A’, and 1” starting at Pj in 2”. These trajectories may be followed 
lor increasing ¢, either for all ¢ > 0, or until some value of ¢t for which they 
pproach the boundaries of their respective regions when the transforma- 

n from Z' to #” becomes meaningless. 

Since the initial values are bounded, db/dt and b(dd/dt) are not greater 
han O(u), so that if P, is not near to the boundary the trajectories will 
remain within their respective regions for an interval of time 7 which is 
Olu). Now by (32) 

h(D)x = K,(b, 4, t) W{<)K(b, 40), 


c 


Fil os » 

h(D)Sx—K(b, dt (= sd ue \K(b,¢,t) = Ke, (35) 
| \ot dt} | ‘i 

where A, is a sum of a finite number of terms containing partial differential 


eiicients of A’ multiplied by differential coefficients of 6 and ¢ with 
respect to t, which will be O(n). 


Let A: be a root of Ai 2) 0 of multiplicity q;- 


J 


Then (35) is equivalent 
tl 


le integral equation 


r— K(t) SYS Sle eA +-C' 
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, 


in which the C),, are the constants of the partial fraction expansion 


aj 
qd. 1 7 
. S Ca 
h(z) a £ (2-2)! 


and the C,; are constants chosen so that the initial conditions at P/ are 


satisfied. Thus C,; and C); are O(1). The summation > is over the distinct 


} 
roots of h(z) = 0. Let the real part of A; be ;, then by III 


Xe < —1/T. 
Any term of the first type in (36) therefore satisfies 
Y QpAj;t| . Y dp—t|T 
C5 tM Cyjlte™. 


If g = 0 such a term steadily decreases, and is reduced to O() in a time 
which is O{T log(1/)}, well within 7,. This may be regarded as the normal 
vase, in which f(z) has no multiple zeros. In the exceptional case when 
h(z) has multiple zeros, q will have non-zero values, certainly not exceeding 
(n*—2), in some terms of (36). Such a term is initially zero, rises toa 
maximum value (omitting the C,;) of (q7')“e~4 at t = q7’, then decreases, 
tending to zero ast + 0. The maximum value, regarded as a function of 4, 
is 1 at q = 0, decreases to a minimum e~"” at q = 1/7, and thereafter 
increases without limit with further increase in g. Clearly some restriction 
on q or T' is necessary if such terms are to be reduced to O(,) in a time of 
O{T log(1/)}. Such a condition would be q; 7’ < C for each q;, with Ca 
constant which is O(1). That is, the time constant of normal modes of the 
linear equation associated with multiple zeros of the characteristic equation 
must not exceed a limit which is smaller the higher the multiplicity of the 
zero. 

If a condition of this sort is satisfied, all terms of the first sort in (36) will 
be reduced to O() in a time which is Of7 log(1/)}, and since there is only 
a finite number of them independent of 1, the same will be true of their sum 

Considering now terms of the second type, we shall have | A,(¢—7)| <m 
where m, is a constant of O(u), at least for a time 7,. For t < 7, we hav 


1 
t t 
P -q—-1 ’ q—-1 
K,(t—7r) erit dr| - | my, . ell dr 
i : (q—1)! / (q—1)! 
0 0 
a ee _\ 
m, 7 | . eit q(T 
(q—1)! J \7 1 
0 
m, T%, 


which is O(4) so long as m, is so. Again, since there is only a fixed number 
of such terms, their sum will be at most O(y). 
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We have, therefore, the result that for any initial point Pj not near to 
the boundary of #”, the difference x—K is reduced to O(yu) in a time 
OT log(1/u)}. In a further interval of the same order S(xv)—S(K) will be 
(Ou), therefore db/dt and b (d¢d/dt), and therefore also Kg, will be O(u?), i.e. 

_may be taken O(u*), and the difference x—K will be reduced to O(u?) 
nd will thereafter remain not greater than O(,"). 

The values of db/dt and 6 (dd/dt) will then be given within O(y3) by the 


utonomous system 


f h i) = 
_ pdb \ cos AU cos(¢- )} (37) 
dd . ‘ y P 

b - ob LA sina b? — ; sin (d t I, (38) 


obtained from (33), (34) by omitting terms which are O(u3). Consequently 
ina further interval of O(1/) the projection of I” on the (6,4) plane will 
not differ by more than O(u?) from the trajectory ['{ of the autonomous 
system which starts from the same point at the beginning of this interval. 

Three cases then arise: 

i) If '} approaches a stable singular point, x approaches the periodic 
solution, the amplitude and phase of which are given in the first approxi- 
mation by the b and ¢ of the singular point. 

ii) If '; winds on to a stable limit cycle, then at the end of the interval 
the projection of the representative point of I", which may deviate by 
Ou?) from the corresponding point of Ij, will coincide with a point of 

nother trajectory [3 of the autonomous system. We may in this way 


btain a suecession of ares Tj, 3, T3,... of the autonomous system, such 


that the projection of [” follows each in turn within O(”) for a time O(1 /). 
Clearly, unless the initial point was near to a separatrix, all these arcs wind 
n to the same limit cycle, and ultimately the projection of [” will remain 
within O(u?) of this limit cycle. 

The solution in this case is an almost periodic oscillation having the 
character of a combination oscillation. The period of the cycle of variation 
tamplitude and phase will be O(1/?) and it will be given with a fractional 
error O(u) by the period of the limit cycle of the autonomous system. 

ii) The trajectory ['} may pass outside the circle b < b,. In this case 
igain the projection of I” will follow a se quence of ares of trajectories of 
the autonomous system, which, if the starting-point is not near to a separa- 
‘nix, will all pass outside the circle. The solution in this case is an oscillation 
‘increasing amplitude which ultimately passes outside the validity of the 


equations of the first approximation. 
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8. Differential equation of the second order 

In the foregoing section it is assumed implicitly that n* > 2. Ifn* = 2 
then A(D) in (4) is a constant which may be normalized to 1. In this 
case 2 = u and S(x) = S(K) throughout. Then (33), (34) are obtained 
without the terms in (S(a)— S(K)), and form a system in b, ¢ only, which is 
altogether equivalent to the original differential equation. This case js 
discussed in (1). 


9. Significance of the order of the differential equation 

The argument in section 7 depends on the fact that terms of O(,:2) and 
higher order are small in comparison with terms of O(ju), and the latter 
small in comparison with unity if » is sufficiently small. At several points 
it is assumed that the sum of a finite number of terms is O() if the number 
of terms is independent of » and each term is not greater than O(y). The 
number of terms in such cases increases with n*, the order of the differentia 
equation, and clearly if n* is large the restriction on p is more stringent than 
if n* is small, in order to ensure that such sums remain small in comparison 
with unity. 

It appears then that small non-linear terms may be more significant in 





equations of higher order, and one might formulate the rough criterion that 
n* must remain a first-order small quantity. 


10. The variational equations } 

The variational equations (37), (38) fall under four main types, according 
as to whether 5 is +1 or —1, i.e. € positive or negative, and whethe 
cosA is positive or negative. There are two further subsidiary cases if cos) 
is so small that the term A cosA 6? is of higher order than the other terms 
retained in the variational equations. 

With suitable changes of scale and origin of ¢ these may be reduced t 
the forms listed below: 

l. e > 0 and cosA > 0, 

b = b(1—b?)—F cos, 
bd F sin d—b(o-+-vb?). 

This is the form obtained under VIII corresponding to (19), (20). I 
includes the symmetrical (v = 0) and unsymmetrical (v > 0) van der Pol 
equation, (1), and the usual case of the phase shift oscillator (v < 0), (2 

2. € < 0 and cosA > 0, 


6 = —b(1+b?2)— F eosd, 


bd F sin d—b(o+vb?). 
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This case is typified by a slightly damped oscillator with non-linear 


restoring force. 
9 ¢€>0 and cosa 0, 
b = b(1+6?)—F cos 4, 
bd F sin d—b(o-+-vb*). 
4,€< 0 and cosA < 0, 
b b(1—b?)— F cos, 
bd F sin 6—b(o-+-vb?). 


The trajectories of 3 and 4 are essentially the same as 2 and | respectively, 
but described in the opposite sense as ¢ increases. 


The subsidiary cases are 
v. ¢ v b b—F cosd, 

bd F sin d—b(o+vb?). 
0. € v b -b— F cos 4d, 

bd F sind—b(o+ vb?). 


Of these 5 is essentially the same as 3, and 6 is essentially the same as 2. 
Case 1 is discussed in some detail in (1) and cases 2—5 are discussed in 
uutline in (2), the nature of the singular points and form of the resonance 
curves being indicated. 

In cases 1, 2, and 6 the equations of the first approximation determine 
the final state of the system for sufficiently small values of ». In the other 
three cases trajectories run outwards for large values of b and from suitable 
initial conditions an oscillation will grow beyond the limits of validity of 


the equations of the first approximation. 


11. Symbolic derivation of the variational equations 

The periodic solutions are determined by the amplitude equation. If we 
suppose that an existing periodic solution is slightly disturbed, so that the 
subsequent motion may be regarded as an oscillation of slowly varying 
amplitude and phase, b and b¢@ being O(y?), then we may neglect higher 
differential coefficients altogether, and take account of the variation of b 
and ¢ in the first-order terms only. The amplitude equation (15) may then 
de written 


G(D)be g(t)(A’ +-1A") 7b Eyre tet (39) 


Within O(u4) with G(D) f(D)+e,9(D). 
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The first term may be written 
elG(D+iw)beld = ci{G(iw) + G' (iw) Dei? 

G(iw)beiot+9) + t (iw)(b+ ibe iot+¢) 

€ Ae ia(§ — ic)be (aut +h)___ Ag inh | ibd etot+d), 
Substitute this expression into (39), putting je} = pu”, separate real and 
imaginary parts after dividing by p?Ae’*, and the variational equations 
(37), (38) are obtained immediately. 


12. General remarks on the method 
(i) The approximate linear equation has been taken as 

{(D—e)? + Q7}h(D)x = 0. (40) 
Neither this form, nor the normalization to Q = 1 is essential. If the 
approximate linear equation is 

{(D?+- O?)h(D) +e, k(D)\x = 0 (41) 

the same procedure may be used, as is done in (2) with Q, normalized to 1, 
If we take z = iQ, as an approximate root of 

(22+ OF)h(z)+€, k(z) = 0 


we obtain for the true value of the root 


€, k(iQ : 
— ~ sae (<1). 
If this is written z= 10,+¢,(M+i1M’) 
within O(¢?) we shall have e=n¢k 
and Q= O,+6, eM’. 
With €,0, = —— and eo = _—> 
Q, OQ 
we obtain eo ala -). 


Thus with (41) instead of (40) the ratio e/e, = M enters the equations as 
a scale factor, and a change of origin of o, to —M’/Q is required to bring the 
variational equations to the standard form. These equations are also 
obtained by the symbolic procedure of the previous section. 

(ii) The relative importance of the non-linear terms has been related, in 
the natural way, to the position of the term in the Taylor expansion of y 
in powers of x. A different assignment of the relative importance of these 
terms could have been made, and the expansions carried out in powers of # 
instead of in powers of B or b. The procedure would have been essentially 
the same, though the results would have lost some of their symmetry. 
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iii) The case e« 0 has not been mentioned explicitly. The non-resonant 
olution fails completely in this case but the resonant solution may be used. 

iv) Fundamental resonance has been considered, but the method may 
be adapted to other cases; e.g. if the forcing term were 2B cos 2wt, with 

| = eo, and x 2b cos(wt+-¢), equations for subharmonic resonance 
f order 2 could have been derived. 

vy) An upper limit to the error made by terminating the series at any 
particular term may be obtained from the remainder after the same number 
fterms of the majorante series. 

vi) In the case when the degree of g(z) is higher than that of f (z)-+-c¢, g(z), 
the problem becomes a singular perturbation problem. The solution fails 
because {(z) tends to zero at infinity and it is no longer possible to choose V 
ess than C(riw)| for all r. If, however, |{(riw)| remains greater than C, for 
r,, say, then with V = C, we may develop the solution as far as the 


4 
tems containing up?” where 2n <1,. If b and ¢ satisfy the amplitude 
equation as far as this term, a solution will be obtained which satisfies the 
differential equation with a residual less than NV times the remainder after 
this number of terms of the majorante series, which may be very small if 


r, is large. 
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RESPONSE FUNCTIONS OF LINEAR SYSTEMS WITH 
CONSTANT COEFFICIENTS HAVING ONE 
DEGREE OF FREEDOM*+ 
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SUMMARY 


The transient response is considered for systems satisfying linear differential 
equations with constant coefficients. Criteria for optimizing the response are giver 
in terms of the response functions 


x 


- 9 


x 
. * (de\? 

L e?dr and L, | ( dr, 

; J \dr 

0 0 

where e¢ is the error at time 7. Expressions are obtained for both Z and L, in terms 
of (a) the coefficients of the characteristic equation, and (6) the frequency respons 
of the system. It is shown how the response due to (i) a step function disturbance, 
(ii) an initial impulse, and (iii) a constant velocity input can be simply related to th 
response in the free motion. 


1. Response functions and optimum response 

WHEN considering the performance of a system (e.g. a dynamical system 
or servomechanism) we are interested in the accuracy with which the output 
of the system follows the input. More precise elaboration of this general 
statement depends upon the particular application, the order of accuracy 
and the sensitivity of the system varying greatly with different applications. 
The systems considered may be mechanical, electrical, hydraulic, or aero- 
dynamic. The terms used here (e.g. forces, equations of motion) are mainly 
based on mechanical systems, but the analysis is, of course, perfectly 
general. 

We shall be concerned, in general, with stable systems, i.e. systems 
which when subjected to any disturbance acting for a finite time ultimately 
return to their initial state. The performance of such a system is intuitivel) 
measured by such factors as its overshoot when subjected to a step dis 
turbance, the oscillatory nature and damping in that case (the transient 
motion), and also by the amplitude of the motion and the resonance peak 
and frequency when the system is subjected to a steady sinusoidal dis- 
turbing force (the frequency response of the motion). Often the designer 
has freedom to choose the value of many of the parameters of the system 

+ This paper forms part of a thesis for which the Ph.D. degree of the University 0! 
Glasgow was awarded. 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 3 (1957)) 
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eg, degree of damping, spring stiffness). The problem of optimization is 
the selection of such values of these variables that the response of the 
system is ‘the best’ or, more often, the most satisfactory for the particular 
pplication. The optimum values will depend on the particular input 
listurbance. Thus at the start we are confronted with the choice of basing 
wr analysis either on the transient behaviour of the system or on its 
frequency response. The nature of this fundamental choice will depend 
on which is the most likely type of input the system will have to deal with. 
Fortunately many systems having a satisfactory transient also have a 
satisfactory frequency response. This is not surprising since the response 
fasystem to any disturbance is governed by the differential equation of 
motion of the system. In fact the transient and frequency response can be 
correlated by a Fourier transformation. Empirical relations are often used, 
based, for example, on the relationship between the peak overshoot and 
the resonance peak, or between the resonant frequency and the transient 
scillatory peak, or between the resonant frequency and the speed of 
response (see (3)). 

We shall endeavour to give simple mathematical criteria for optimizing 
the transient response. Many such criteria have been used in recent years, 

inly for systems subjected to input step disturbanees. In references (2) 
nd (3) the intecrals 

I, |}edr and J, 


0 0 





where ¢ is the error) are minimized. These criteria break down when the 
response is oscillatory, since an overshoot decreases the value of J, and J. 
In (4) and (5) the criterion for optimizing the transient response is that 


should be a minimum. This criterion is widely used and can be readily 
handled, either analytically or by a computing machine. However, in 
yeneral it leads to a slightly overdamped response often with a large 
indesirable overshoot. To overcome this Mack considers (6) the integral 
I, | re? dr. 
0 
the system making J, a minimum gives a very satisfactory performance, 
ut the formula for J, is often troublesome to evaluate. 


In reference (7) 
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the integral of time-multiplied absolute-value of error (ITAE) is minimized, 
The ITAE criterion works very well for response of zero-displacement- 
error systems to a step input, but not so well for zero-velocity-error systems, 
It may indeed be thought to be over selective in its choice of optimum, 
distinguishing too sharply between the optimum and systems near the 
optimum. From the form of the integral it allows sizeable errors for smal 
values of 7; this is seen especially in its choice of optimum response for 
zero-velocity-error and zero-acceleration-error systems. While J; is easy 
to calculate, it is difficult to handle analytically, and thus it is hard (if not 
impossible) to extend results of the ITAE criterion to high order differential 
equations. 
We shall consider here the response functions 
rs v3 
L=|€@dr ad L,= | E 


dr 


9 


dr. 





. 


0 0 
As stated above a system based on minimizing L often has an undesirable 
overshoot. Often Z has a rather flat minimum as the parameters of the 
system are varied, the value of this minimum being relatively insensitive 
to small changes in some of the parameters. The system having the most 
satisfactory performance is impossible to define precisely in mathematical 
terms. For a system with a large overshoot de/dr will be large at some 
time and thus L, will be correspondingly large. Our aim will therefore be 
to choose values of the parameters so that LZ is near its minimum value 
while LZ, is considerably lower than its value at L 


min* [hese systems will 


be of a less oscillatory nature than those given by the mean square optimum. 


2. Derivation of simple formulae for the response functions L and |, 


We consider a system for which the equation of motion is 


a, D°x+a,_, D"x+a,_, D"—*x+...+a, Dxr+a x = f(r), (I 
d 
where D=-—, 
dr 
Big, By) Bey. a, are constants (a, + 0) and f(z) is some known function. 


The initial conditions are given by 
t=@, De= Dx, Pa= Dr ch, D2 = Ds, attr =t 
The characteristic equation corresponding to (1) is 

F(A) = a, A"+a,,_,A" +4, _.A"-?+.... tay = 0. (2 


n-2 
We shall be concerned only with stable systems for which all the roots 
of (2) are negative or have negative real parts. 
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If there were no lag in the system x would equal (1/a,))f(7) at all times. 


The error e is given by 


L f(r)—x. (3) 


a 


é input—output 


As stated above we shall derive formulae for 


: oe a we 
L = a dr = | & fia)| dr (4) 
{ Pi 7) ; 
0 0 
rTde]?,  [fdx 1 : 
‘ai 8 dr = — “(fr Ir. 5 
id Ly | Z| | E a! ( ] ; ) 





0 0 
We shall now consider various forms of the functions f(z) and the corre- 


sponding formulae for the response functions LZ and J,. 


3, Free motion: f(7) = 0 
: rs 
We write L = | 2dr L, = | [D*x}? dr. (6) 
0 0 
Putting f(r) = 0 and multiplying (1) successively by x, Dz,..., D"-1x” and 
integrating from 0 to oo we obtain simultaneous linear equations for 
[,,..-, L,-1 in terms of 25, Dxp,..., D"-125. We find 
L l B, Ge @ Ge... (7) 
iy T-1 B @& Gs @& « « « 
. me wm « 4% 
6B, 0 a a, 
| 
nd a I Bo Gy Gy @ « ss (8) 
Rid Bs Gy Gy ag 
B, @, Gy @, 
8. Ge Gy 


where 7, _, is the (n—1)th test function (see (8) and (9)) given by 


4 i | ( 
- |Q@, a3 as G@, . (9) 
Ay Ap Ay Ag 
U0 @ a3 a; 
O Ay Gy Ww 


he roots | and the determinants in (7) and (8) are of the nth order while the deter- 





inant in (9) is of order (n 


-1). 
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The 6’s are quadratic functions in the initial conditions 2p...., D*-z,, 


Thus 
By = 2% (4, %+a, Dxy+a, D?x,+a, D¥x,+...)— 
— Dx,(4a, Dxy+a, D?x, +a; D®x,+a, D4x +...) + | 
+ D*x,(4a; D?x,+a, D°x,+a, D4x,+a, D®2,+...)— 
— D8x,(4a, D3x)+a, D4x 9 +a, D'x) +449 D®xy+...)+..., (10 
By = $a) 22+ Day(a. Dxy+a, D?xy+a, D8x,+...) 


— D*x,(4a, D?x,+a; D'x,+a, D'x,+...)+ 


+ D3x5(4a, D®x,+-a, D*x)+-a, D'x)+...)—..., (1 
Bs = Da o(a) 2+ 4a, Dxy)+ D?x,(4a, D2x,+a, D®x,+a; D4x,+...)— 
— D®x,(4a; D?x.+-a, D*x)+a,D*®x,+...)4 
+ D4x5(4a,D4x.+-a, D®x +a, D8xy+...)—...; (12)? 
By = —hay Dai + D?24(a9 2%) +a, Dxy+ 4a, D?x)+ 
+ D8x(4a, D3x,+-a; D4x,+a, D®x,+...) 
— D*x,(4a, D4x, +a, D*®x,+a, D®x,+...)- 
+ D5x,(4a, D'x,+-a, D®x,+-A4) D'xy+...)—.... (13 


These equations simplify considerably when the initial values of all except: 
one Dz are zero. 

As shown in (8) and (9), for stability (with a, > 0) all the test functions} 
must be positive. The critical stability criteria are 7, > 0 and a, >! 
When a, vanishes, the characteristic equation (2) has a zero root and th 
system is neutrally stable. When 7),_, vanishes, the system has a pair 0! 
equal and opposite roots. In both of these cases if the system is displacet 
from its original position it would never return to and remain in that 
position; thus L is infinite if either a, or 7,,_, vanish, as shown by (7). 

We see that ZL and L, are both second degree expressions in the initi 
conditions. An alternative expression for Z, can be obtained directly fro1 
L by replacing x) by Dxy, Day by D?xp,..., D”— 25 by Dx, in the expressiol 
for L. 


4. Derivation of integral formulae for the response functions it 
terms of the frequency response spectrum of the system 
As stated above, the transient and the frequency response can be corre 
lated by a Fourier transformation. It is of interest to relate the respons 
functions to the amplitude and phase of the frequency response charac: 





teristics of the system. 
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RESPONSE FUNCTIONS OF LINEAR SYSTEMS 
As shown in (10), if 


G(w) = | e-t72(r) dr, 


(14) 


) 


l 
tnen A ( T) € ‘oT G(w) dr. 


7 
. 
x 


provided that { a(7) dr is absolutely convergent. We shall consider only 


the special case in which a(7) = 0 for 7 < 0. 


In the {ree motion we have 


— ~. A,T 5 
t= > Aw*, (15) 
r=1 
here A. ( | to n) are determined by the initial conditions and A, is a 


tt of the characteristic equation (2). From (14) and (15) we have 


G(w) Ss " A, 


a—t iw A, 


rovided the system is stable. On substituting for A, we find 


( 1,(iw)-+-qo(tw)?+...+-¢ iw)"—1 . 
Bis Yo l é Yo ‘ Gn 1( a ae (16) 
Ag -A,(tw) Ao\(tw)*--... a,,(v@)" 
where 
g A414 +a »D gto tA Dr s—ly, (s 0 to n—1). (17) 
From Parseval’s theorem, for stable systems, 
lf Fo 
G(w)|? dw | x?(7) dz, 
0 0 
since a(7r) 1s real. 
Thus, in the free motion 
lf lg Jo w"--g w4 ... |? +w | ¢ g,w?+... |? 
L \Jo—Ie Ja '+ te +E 4, (18) 
T ) |Qo Ay W aw? coal t w*| a, a,w* Teco] 
u 
where g, is given by (17), s 0 to n—1. 
Equation (18) is an extension of the formula given in (6). For sufficiently 
g : 


large values of w the integrand tends to g,_,/(a, )*. Thus the integral 
converges like 1/w. 
Now from (14) we see that |G(w)| represents the frequency spectrum of 


the given system. If the system is given a periodic disturbance a, e'®’, we 


find from (1) that 


i Eee 


Ay +d, iw—d, w*—a, w+... 


0 
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The amplitude (or dynamic magnification) M and phase advance N ar 
given by 


] . 
= 5 {[@9—@g w? +a, w4—... ?+-w{a,—a,w?+...] (19 
Pe ao 
3 5 
- a, @—G,.w”" 4. wr" —.... 
and N a Big fe esl Bis —} (20 


Ay — A, w? +a, w*— 
0 2 — 


re 


We see that the expressions on the right-hand sides of equations (19) and 
(20) are precisely those that occur in the denominator of (18). Large 
amplitudes will occur when the damping is small and the forcing frequency 
is close to one of the natural frequencies of the system. From (19) it is set 
that if the dynamic magnification is large over a range of frequencies, the 
integrand in (18) would be expected to be correspondingly large and | 
would be well away from its optimum value. 


5. Response of a system with no initial displacement 
We consider a system with initial conditions 
a 
the system satisfying the equation of motion (1). 
(i) Step function disturbance 
f=090 forr <0; fj=f, forr> 09, 
where f, is a constant. Writing 


So 


H r~— 
Ay 
we see from (1) that the response of the system following a step disturbance 
is identical in form with that in the free motion for which 
eo —< gud De, D*x....., D"-1xy, are zero. 
a 
0 


The errors in the two responses are the same at any given time and the 
values of Z and L, can be found from (7) and (8). Thus 


"2 A) 
L- fo Mm G % UW + + +) - 
ee 
=1y T’, 1 Ao ay as as 
0 Ay As Ay 
0 QO a, As 
2 99 
i Jo ay Ay Ag Ag . - «|e - 
4} T 
=Uptn-1| % 43 Gs ay 
My Ay Ay a 











*) = Da, = D*x, = ... = Dx, = 0, (21) 
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it follows that the system having the most satisfactory response to a step 


fnction will have the most satisfactory response in the corresponding free 


motion. 
19)| 


19) and 

Large 
quel C\ 
1S see! 
ies, the 


and | 


rbar 





} 








ii) Response to an initial unit impulse 
Itis readily shown that just as the unit step function is the integral of the 
nit impulse, so the response to a unit step disturbance is the integral of 


the response to a unit impulse. The response functions can thus be obtained 


lirectly from the response functions for a step function disturbance. We 
find L, (unit impulse) = L,.,, (unit step disturbance), (24) 
where L, L. 

iii) Constant ve locity input 

] 0 fort 0; j=h, -fyt for 7 > 0, 

vhere f, and f, are constants. 

When the transient has died away, the steady motion is given by 

2 ar+b, (25) 

where, from (1), f, =a, a and f, = a,a+a,b. (26) 


Writing r’ x—ar—b 


we see from (1) that the response of the given system is identical in form 
with that in a free motion with 2, b, Dxy a, and D*zp...., Dp *-tz, 
zero. In particular when f,/a, = f;/a,,6 = 0 and the equivalent free system 
has Dy f D"—x, zero. We see that such a system 
iltimately has no statie error since x’ > 0 and there is no position error 


a, and 2, D*a 


0 


with constant velocity input (since b 0). From (7), 


] E i Ag A, A, a (27) 
Fi yo 
Wp 1 1 As ay ag a- 
a, Ay As 4 
Ao 0 ay a3 


Similar formulae can be deduced for a constant acceleration input. 


6. Extensions to the above analysis 

\s is well known, the response of a linear system to an arbitrary distur- 
bance can be simply related to its response to unit impulses. When the 
nput is expressed as a Fourier integral the response functions L and L, 
can be found by an extension of the method given above, using Parseval’s 
theorem. The response functions for a finite square wave disturbance could 


be obtained in this manner. 
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7. Response functions for unstable systems 

The response functions defined above apply only to stable systems, |; 
sometimes happens that the system possesses one slight unstable mod; 
(e.g. an aircraft with spiral instability). In this case we are primar] 
interested in the behaviour of the system shortly after the disturbane 
takes place. The response functions can be used to give a measure of the 
response of the system, the upper limits of the integrals for Z and L, being 
taken to be some convenient finite time. 


REFERENCES 
1. G. J. THALER and R. G. Brown, Servomechanism Analysis (McGraw-Hill, 1953 
2. R. C. OLpDENBOURG and H. Sartorius, The Dynamics of Automatic Contr 
(ASME, 1948). 
3. P. T. Nims, ‘Some design criteria for automatic controls’, Trans. AIEE 7 
(1951), 606. 
4. N. Wiener, Extrapolation, Interpolation and Smoothing of Stationary Time Seri 





(Wiley, 1950). 
. A.C. Har. The A nalysis and Synthe Sis of Linear Servomechanisms (MIT, 1948 


1? 


6. C. Mack, ‘The calculation of the optimum parameters for a following systen 
Phil. Mag. (7) 40 (1949), 922. 
D. GrawaM and R. C. Laturop, ‘The synthesis of optimum transient respons 
criteria and standard forms’, Trans. AIEE 72 (1953), 273. 
8. A. Hurwitz, ‘ Ueber die Bedingungen, unter welchen eine Gleichung nur Wiirzel! 
mit negativen reelen Theilen besitzt’, Math. Ann. 46 (1895), 273. 
9. R. A. Frazer and W. J. Duncan, ‘On the criteria for the stability of sm: 
motions’, Proc. Roy. Soc. London, A124 (1929), 642. 
10. H. 8. Carstaw and J. C. JAEGER, Operational Methods in Applied Mathemat 
(Oxford, 1941). 


“I 











CO 


The é 
set up 
describe 
comple: 
determi 
for poir 
descril y 
senting 
basis fe 
polynot 
time w 
given fi 


1. Th 


1.1. 
LET u 


of deg 
quant 
the ec 


where 
By 


obtail 


Th 
some 
To 


=4 


[Qua 


5092 
































AN ELECTROLYTIC TANK AS AN ANALOGUE 
© mo! COMPUTING MACHINE FOR FACTORIZING HIGH 





‘Imari 
CODER D —_?, : 
rbat DEGREE POLYNOMIALS 
e of 
1 bei By 8S. K. £P (Imperial College, London) 
Received 25 June 1954; revise received 15 January 1956] 
SUMMARY 

The analogue representation of a polynomial by means of a potential distribution 
ll. 195 set up by point electrodes in a uniform sheet of conducting material was first 
Cor lescribed by Lucas. The electrodes are suitably placed along the real axis of the 

mplex z-plane after a region in which the roots may be expected to lie has been 
IEE 7%\ determined. The roots can then be found by exploring the electrolytic tank surface 


r points of zero potential gradient. In addition to the double-layer type which is 


Se lescribed fully by Cherry et al., a single-layer transformed electrolytic tank, repre- 
ting the infinite z-plane, was employed throughout the problems explored as a 
, 1943). basis for this investigation. By means of simple measurements, all the roots of a 


syst polynomial as high as the 16th degree could be located in between 2 to 4 hours’ 
with an overall accuracy of about 4 per cent. For special cases a method is 


given for setting up the potential field by electrodes placed along the imaginary axis. 


Wit 1. The physical representation of algebraic functions 
1.1. The algebraic function analogy as described by Lucas 
Let us consider the equation 


f(z) a,2z"+a gn 


1 » - 
a +-...+a,2+a,) = 0 (1) 


f 


[degree n with real coefficients and let us take n+-1 arbitrary unequal 
quantities x,, 2,..., ®,,,, each determining a point Z on the real z-axis of 
the complex z-plane. Let us now form the auxiliary polynomial 


n+1 
F(z) Il (z—a,) (2) 

s=1 

where F(z) is of higher order than f(z). 
By dividing f(z) by F(z), where the coefficients of each are real, we 


tain a rational funetion G(z) which can be expanded in partial fractions: 
G(z) ——. (3) 


rhe residues k,, ky, ete., can easily be determined. They are all real, 


some positive and some negative. 





lo determine the positions of the roots of f(z), the complex plane of 
«ty is represented on a uniformly conducting sheet (electrolyte), 

Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 3 (1957)] 
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and each of the partial fractions is represented by feeding a current of 
strength k, into the sheet at the positions 2,. A system of current flow js 
then set up which characterizes the function G(z). 


1.2. Conditions at a root point in a potential field 
Let us consider a point P(z) in a potential distribution 
; n+1 
V (z) > k,log(z—€,) (4 
s=1 
due to a number of line sources of strength /, at positions €,. The conditio 


that V(z) be ‘stationary’ near a point z, is dV /dz, = 0, that is, 


Sn 

— (2 &.) 
which can be seen from equation (3) to be identical with G(z) 0. This 
same condition constitutes the solution of f(z) 0 in (1). 

Then if G(z) be thus represented by means of the gradient of a potentia 
distribution, the positions of the points z, will be the roots of the propose 
equation f(z) = 0. 

The positions of the points z, in the potential field can be traced out by 
means of a single probe. This makes use of the fact that at the point : 
where dV /dz = 0,i.e. the point with zero potential gradient, there is a saddl 
point of potential, as follows from Dirichlet’s theorem. 

To investigate the behaviour of the equipotentials in the vicinity of z 
let us put 


so that V(z) V (z)+Az). 
Then expanding (5) in a Taylor series, we have 

V (zy +Az) = V(z9)+AzV'(z)+4(Az)?V"(z,) + O(Az)3. 
Hence the equipotentials in the vicinity of z, are given approximately | 

(Az)? = constant, 
that is, either the hyperbolas 
(Ax)?— (Ay)? constant 

or AvAy = constant. 

Thus equipotential lines around z, occur in hyperbolic form. The limiting 
equipotentials (the asymptotes of these sets of hyperbolas) are straighit 
lines that intersect at z,. Moreover, since z, is a saddle point, A B and (1 
in Fig. 1 are loci of potential minima and maxima respectively, intersecting 
at the root point. 

On the electrolytic tank, current-carrying electrodes of small diamett 


set up a distribution of potential which is peculiar to the given polynomi 
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The potential maxima and minima are traced out; their intersections will 


ie the positions of the root points of the polynomial. 











Equipotentials 


Fic. 1. Seetior iew of a saddle point of potential 


2. The electrolytic tank method 


It is often extremely difficult to factorize polynomials exactly, but there 


re approximate methods for numerical solutions (2)—(8). 
The following is an electrical method based on the practical experience 
the electrolytic tank analogue found by Boothroyd, Cherry, and Makar 
9), and is chnique of the greatest value for representing a complex 
It is applicable, in many cases, for solving both the real and complex 
toots Ol rational algebraic equations of any order that have real coefficients. 
The equipment « mployed is simple and inexpensive compared with other 














372 So. B. IP 


mechanical and electrical computing machines (10), (11), (12), and further. 
more is easy to operate. 


2.1. The single-layer electrolytic tank 

By the bilinear transformation, the whole z-plane above the real axis js 
conformally transformed into a unit circle, and is represented by the region 
within the single-layer tank. The relation between z and the new complex 
variable ¢ is defined by 


t(= ——. 


In practical form, a Perspex ring 1 em thick and 40 cm inner diameter 
fixed on to a plate-glass disk, forms the wall of the single-layer electrolyti 
tank. The real axis (—o to «) of the z-plane now becomes a circle of 
radius $ and centre at 37 in the ¢-plane, and forms the boundary of th 
tank. Since the complex roots of a polynomial with real coefficients ocew 
in conjugate pairs we need deal only with one-half of the z-plane (either 
above or below the real axis) and the tank gives us the whole of such 
region. Note that the coordinates which are marked on the plate-glas 
must be curvilinear instead of rectangular. 

Copper sulphate solution (concentration 4 to 5 g to 1 |. of distilled water 
to a depth of 1 em is used as electrolyte. Operating currents, of the order 


of milliamps, are fed from a 500 ¢/s stabilized current source; positive and | 


negative currents can be taken as being in opposite phase to each other. 
In cases where the roots have large numerical values, or all are offset t 
one side of the imaginary axis, we can bring them into the favourable regio 
2 to —2 on the real axis by a suitable linear transformation. There is 1 
need to apply this transformation to the equation itself, but only to th 
roots when they are found. 
For example, suppose the roots lie between « = —150 to —50. Then the 


‘ansformation is 
transformation is 2+100 


so that a root z, in the coordinates of the tank gives 25z’— 100 as a root 0! 
the equation. 

If we put the electrodes at —150, —125, —100, —75, —50, this would 
imply that they are placed at the points —2, —1, 0, 1, 2 in the z’ case. 


2.2. The positioning of the electrodes 

The positioning of the electrodes 2, 2», etc., is arbitrary. But, since th 
electrode field intensity at a point is inversely proportional to the distance 
between the point and the electrodes, the latter should be placed as near 
the root points as possible. If it is possible to find a limiting region S on the 
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‘urther.| real axis, inside which lie all possible roots of a polynomial, this will be the 
ntimum region in which to place the electrodes. 


One method of determining the region S is by plotting f(z) against the 





real a-axis: in the vicinity of the roots, the variation of the function is more 
| axis gradual than it is at a distance from the roots. When a polynomial is made 
 TZION | fee from real roots, as was generally done (by standard methods) so as to 
omplex} reduce the polynomial to the lowest possible degree, the graph will not 
_ | touch the real axis and the region S can be taken as lying between the two 

steep slopes of the graph as shown in Fig. 2 (a). 
aiden \fter the region S is determined, n+-1 real values x,, x, etc., are arbi- 
trolytie| ttatily chosen from it (section 1.1). When possible, it is desirable that 
ircle off these 2+ 1 values should be distributed evenly within the region. Further 
of thet they should be convenient values, otherwise unnecessary complications 
a aie’ nd labour may be entailed in the solution of the partial fractions that 

either} follow. 


When the roots of a polynomial occur alongside the imaginary axis, f(z) 


te-clag | can be formed into a new function H(z) = f(z)f(—z). By plotting H(z) 
gainst the imaginary axis, a similar limiting region S can be obtained 
nthe imaginary axis. In this way potential distribution can be set up by 
ectrodes along the imaginary axis. 


— \ general method for obtaining the boundary values of the roots is to 





other. | Plot f(z) against both the real and imaginary axes; hence the roots can be 
fiset | ited approximately within the rectangle bounded by the four limiting 
y ree ints A, B, C, and D, as shown in Fig. 2 (b). 
satias 2.3. A quick and accurate method of searching for the root points of a poly- 
i nomial on the electrolytic tank 
On the electrolytic tank a single moving probe is used to measure poten- 
aoe tials with respect to an arbitrary point on the tank surface. The probe 
sweeps the tank, firstly in closely-spaced parallel lines and then similarly 
the normal direction, until the whole area of the tank has been 
sia explored. 

The potential at various positions can meanwhile be observed on the 
~— galvanometer. When the potential reaches either a maximum ora minimum, 
cast the position on the tank surface is recorded on a sheet of graph paper 

isually having the same size as the tank surface. At the end of the search 

throughout the tank surface, the locations of the saddle points can 

nce t mmediately be found where two loci of potential maxima and minima 
istat ntersect. 

near t Then within appropriate regions in the vicinity of the intersecting points 

equipotential lines are carefully traced out. The saddle points of the equi- 
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Fic, 2. Approximating the locations of the complex roots of a polynomial by 
an area bounded by four limiting points 
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notential lines in cartesian coordinates will give the complex values of the 
ots of the polynomial represented by that potential distribution. In 
ractice, the position of a saddle-point cannot be located precisely by means 
fpotential field configuration; however, a null region is found and the exact 
sition of the saddle point can thus be taken at the intersection of the 
symptotes of the equipotentials as shown in Fig. 1. 

(a) Care must be taken when searching for the potential saddle points. For 

example, Fig. 1 illustrates a sectional view of a saddle point of potential. 





fin this case the probe travels along Mz, N and Pz, Q, that is, along an 
ipotential, neither a maximum nor a minimum will be recorded. It is 
sometimes necessary to make further explorations in a direction of 45° from 
ch of the original paths of probe movement, that is, along Az, B and 


Cz,D. 





/ 3. Practical examples 
| The two examples given below were solved by the transformed single- 
wer electrolytic tank, on which electrodes were positioned along the real 
Example 1) and imaginary (Example 2) axes. 
Example | is due to Hitchcock (7) and was also solved by Lin (8). The 
‘ts are well distributed, which is one of the most favourable conditions 
solution by means of the electrolytic tank. A close degree of accuracy 
vas obtained. 

By placing the current-carrying electrodes along the imaginary axis, 
Example 2 was solved with the same procedure given in Section 2, in spite 
‘b ) {the high order of the function. 

Example | 

Here 


3-01 227 + 3-22526+ 1-021225+ 6-9862z4— 21-8892? 


8-11022--5-9012z-+ 23-889: (8) 


is plotted against the real axis in Fig. 3 from which the electrode region 





Sis chosen between — 2-0 and 2-0 with the x’s at 2-0, 1-5, 1-0, 0-5, 0, —0-5, 
1-0, —1°5 2-0. The values of the residues are found to be 
k, = 0-726012 k, ~9-468622 k, = 10-877688 
hs 1591714 I 10-613733 ky —12-626996 
ky = 4485866 Ke 9-310044 ky = 7:280400 
Li 
The values of the electrode currents are given in Table 1. The loci of 


potential maxima and minima, and the locations of the root points, are 
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shown in Fig. 4, whereas Fig. 5 is a corresponding drawing in which , 
physical picture of the behaviour of the equipotentials is traced out for the 
purpose of illustration. 


F (2) 
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L : 
Fic. 3. The graph of 
f(z) = 28—3-01227 +4 3-2252z°+ 1-021225+ 6-986224— 


— 21-8892° + 8-1102?+ 5-9012z-+ 23-889 


TABLE l 


Electrode Xy Xe Xs Xs Xs Xg Xq Xg Xy 
Position 2°0 1°5 IO 0°5 ° 05 1°0 1°5 2°0 
jsource I 6:2 14°6 15'0 10°0 


Current (mA); . , : 
\sink 218 13°05 | 12°8 | 17°4 





A comparison between the measured results and those given in (8) is 
shown in Table 2. 


TABLE 2 








Values given in ref. (8) 0°48-+ 10°804 1'043 


I 
Measured values I° 


*515-+10°618 r51-+21°5 
512 


5 ' 
10-610 1°504+71°566 | —o-483+10°79 | —1:047LH'O 
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The time required for locating the roots by the tank method is approxi- 
itely 3 hours and the overall accuracy is about 1 per cent. 
It should be pointed out that in Fig. 4 the loci of potential maxima and 


ma do not agree accurately with the configuration of the equipotentials 
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f Example | on a single-layer electrolytic tank 


the vicinity of a root point. The explanation is that the tank surface 


searched by a rapid but approximate method by the searching 


robe. The probe was guided along the orthogonal lines on the graph paper 
tthe bottom of the tank surface usually at 0-5 em intervals. The maximum 
nd minimum potentials, which are in fact the points of tangency to poten- 
il lines shown in Fig. 5, were observed and their positions recorded. 
When a gradual curvature occurred along the path of measurement a longer 
wk was recorded than for a sharp curvature. The loci of potential 
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maxima and minima are lines constructed by joining these long or short 
marks, with correspondingly less or greater accuracy. 
It cannot be expected, therefore, that these loci coincide very accurately 


with the actual configuration of the equipotentials. Also it must be noted 





Fic. 5. Potential field representation of Example 1 corresponding to 
Fig. 4; lines are equipotentials 


that the figure is produced directly from the experimental chart without 
modifications. Although these loci merely gave the clue to the regio 
in which the final potential measurements should be made, the inte! 
sections of these lgci do, nevertheless, give the root points accurately. 11 
the next example attempts were made to locate the root points directly )) 
the intersections of the loci without final potential measurements being 
made, 
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Example 2 
Here 
f(z) = 28+ 2-40215+ 9-6392"4 + 16-922213 4 34-673212+ 45-46521 + 
+ 60-264219+ 59-01829 + 53-97328 + 38-64327 + 24-37726+ 
+ 12-1492°+- 5-01724+- 1-57323 + 0-3682?+ 0-056z+0-00427, (9 


Fig. 6 shows the graphs of f(z) and H(z) = f(z) f(—z) plotted against the 
x- and y-axes respectively. From this it is seen that the roots will occur 
within the regionz = 0, —0-5,andiy = 0, 1-6, that is, in the neighbourhood 
along the y-axis. In order to set up current flow by electrodes along the 
imaginary axis, H(z) is dealt with instead of f(z); in this case, 2n+ 1 electrodes 
are to be used, including one at z = 0. The values of the corresponding 
currents can be obtained after determining the partial fractions: 

: H(z) a _*y : _ he es (10) 
(z—ty,)(z—tyo)...(Z—tYona,)  (2—ty,) | (z—tye) 

Now, because the 2” roots of H(z) are symmetrical about both the z- 
and y-axes, only one quadrant of the complex z-plane need be used for 
locating the roots, and the number of electrodes can also be reduced to 
n-+-1; the value of the current at z = 0 is one-half of that originally 
obtained from (10). 


Electrode positions were chosen at y = 0, +7(0-1, 0-2, 0°3,..., 1-5, 1-6). 
Only the loci of potential maxima and minima are shown in Fig. 7, and 
the positions of the roots were determined by the intersecting points of 
the two loci. The theoretical and measured values of the roots are compared 
in Table 3. 


TABLE 3 














Theoretical 0°210-+10°104 0°185-+10°315 0°174-+10°490 0°162+10°674 
Measured —or1g1+10°094 | —0'175+10°310 | —o-149+10-478 | —o-148+10°673 
Theoretical 0°148+-10°917 o-120+21°14 o-105-+71°31 07098 + 21°51 
Measured —O'IlI+10°917 o116+71'11 —0°095-+11°32 0°078+11°50 


cates | 1 





The time required to locate all the eight conjugate pairs of roots was 
2 hours with an accuracy of about 5 per cent. 


4. Accuracy and advantages 

The accuracy of solutions obtained by this method depends to some extent 
on the experience and judgement of the operator as well as on the intrinsi¢ 
accuracy of the apparatus. However, the overall accuracy of the values of 
the roots obtained throughout the investigation was about 4 per cent. This 
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isconsidered to be a satisfactory result, since the root values obtained may 
be regarded as first approximations which may then be improved to any 
lesired degree of accuracy by familiar numerical methods. 


7-01 mA 





Fic. 7. Root locations of Example 2 on a single-layer electrolytic tank 


Three main advantages result from the method described here, particu- 
larly when polynomials of high degree are involved: 

|) Provided the electrode arrangement is suitably selected, roots of 
polynomials of any degree can, in principle, be located accurately with 
approximately the same speed. 
2) All roots can be clearly visualized at once; this avoids the possibilities 
ofroot divergencies which are common if a solution is performed by ordinary 
humerical methods. 
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(3) When electric network problems are dealt with by means of the tank 
(9), pole and zero positions can be found directly on the tank surface by the 
present method. 

When a polynomial is dealt with by means of the tank, the major part 
of the time is spent in the preliminary evaluation of the function and the 
subsequent determination of the residue values for the electrode arrange- 
ment by substitution in the relation given by equation (3). In the actual 
measurement, approximately half an hour was required to locate one root 
on the tank when equipotentials in the vicinity of the roots were traced out. 
The tracing of equipotentials may not, however, be necessary when only a 
first approximation is required. For then the locations of the roots can be 
determined with sufficient accuracy from the intersections of the loci of 
potential maxima and minima, as in Example 2. The nearness of the root 
points obtained in this way to those obtained by the tracing of the equi- 
potentials is evident from Fig. 4. It may be remarked that the time required 
to find all the roots from the intersections of the loci of potential maxima 
and minima is no more than | hour. 

The accuracy of the results described is limited by margins of experi- 
mental error imposed by the equipment used. This error arises in deter- 
mining, usually at the intersection of the asymptotes of equipotentials, the 
exact position of a saddle point. 

Four possible sources of such error are as follows: 

(1) The degree of accuracy in the measurement of the current values 
required to set up the potential distribution was dependent upon the 
sensitivity and calibration of the milliammeter, variable resistances, power 
supply, ete. Although great care was taken, fluctuations of current values 
due to polarization at electrodes were inevitable in a long experiment. 

(2) The degree of accuracy of measurements was dependent upon the 
sensitivity and calibration of the valve-voltmeter. 

(3) In the actual determination of the equipotentials, the use of 0-5 mm 
diameter platinum wire as measuring probe limited the number of potential 
measurements (for tracing equipotentials) that could be made within a 
unit area around the zero, and consequently the certainty of definition of 
the null region as provided by the tracing of the equipotentials. 

(4) The thickness of the probe may also have affected the degree of 
accuracy with which, within the null region as defined by the tracing of the 
equipotentials, the exact position of the root point was determined. 

In addition, the equations dealt with in the foregoing experiments include 
many of a nature that gave rise to unfavourable conditions for solution by 
the method described. Such unfavourable conditions include, for example, 
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when the roots are (a) very far from both the 2- and y-axes; (b) cluster 
together;+ (c) vastly different in moduli. 

In these experiments, despite the use for each equation of a number of 
lectrode positionings in order to find for each the optimum electrode 
positions, the residues obtained exhibited vast (highly unfavourable) ratios. 

In the next section it will be shown that the greater the ratios exhibited 
y the residues, the more unfavourable the conditions for the method 


lescribed 
5, Justification of the choice of arbitrary roots x,, x, etc., indicated 
by the values of residues 


The determination of the values of the residues from the partial fractions 


btained by expanding G(z), as in equation (3), gives the values for the 
points x, giving rise to the potential distribution corresponding to the 
polynomial f(z) of (1). The process is straightforward yet sometimes 
borious, particularly for high degree polynomials. As the roots in equation 
2) are chosen arbitrarily, many different electrode arrangements can be 
btained, but only that close to the ideal condition will be of interest. The 
leal condition is that the electrodes should be close to the root points: all 


| 


the electrodes should be placed within the region S chosen, and the ratio 
the distance between the furthest root and the axis where the electrodes 
re positioned and the distance between electrodes should preferably be not 
ore than 4 to 1. Because the ratio does not exceed 4:1, the ratio of the 
ighest to the lowest value of the residues will generally be of the order of 


t more than 25 to 1. It is found that the lower the ratio between the 
ilues of the residues, the more favourable the conditions for obtaining 
itions. When a set of residues is obtained, one can then tell whether such 
rangement of electrodes is suitable. If not, fresh electrode positions are 
hosen. In this way, electrode arrangements can be justified by observing 
the values of the residues until a satisfactory set is obtained. This eliminates 


i lot of unnecessary trial on the electrolytic tank. 


When a set of residues is obtained having values of vast ratio, two possi- 
ilities may exist. It may be either that electrode positions were not suitably 
hosen or that the roots are unfavourably distributed. There is no rule to 

prove the positioning of the electrodes, but experience of judging the 
residue values is a great advantage and this may help to determine whether 
such a polynomial can be solved by the tank method at all without the tank 
eeding to be set up 

W gether it is impossible in practice to differentiate one 


1use tl rrent intensity vanishes in their vicinity and only the 


ed. 
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6. Conclusion 
The application of this method for factorizing polynomials of higher order 
than the 16th, which is the highest attempted in this work, is not impossible; 


and theoretically the method can be extended to solve equations of any 


degree. It is hoped that the electrolytic tank will prove extremely valuable 
in the analysis and synthesis of electric networks and theoretical analysis 
of those physical and engineering problems which may be represented as 
polynomials. 
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